{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. 标准线性回归(LR)\n",
    "$$\\mathbf{y} = \\mathbf{X}w + b$$\n",
    "$$平方误差：\\sum_{i=1}^m (y_i - x_i^Tw)^2$$\n",
    "$$导数为零求得：\\hat{w} = (\\mathbf{X}^T\\mathbf{X})^{-1}\\mathbf{X}^Ty$$\n",
    "目标：__最小平方误差__。  \n",
    "注意，需要先判断矩阵是否可逆，判断方法：行列式不为零。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-12-11T12:14:48.106709Z",
     "start_time": "2018-12-11T12:14:47.275313Z"
    },
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAD8CAYAAABXe05zAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl4VOXZP/Dvk8mQTEAyIGEb9jWyiJFUQLQWXKAumLq8dWu1tXXpqvKj4isqWH8FpWpta/uWt7WLWAsKRnABRXCjYg0mQAIJCrINWxQmLJkkk5nn/SM5YTJzzpxzJnPOLHw/18V1JZMzZ55D4J5n7nM/9yOklCAiovSRlewBEBGROQzcRERphoGbiCjNMHATEaUZBm4iojTDwE1ElGYYuImI0gwDNxFRmmHgJiJKM9lWnLRHjx5y0KBBVpyaiCgjbdy48UspZYGRYy0J3IMGDUJZWZkVpyYiykhCiN1Gj2WqhIgozTBwExGlGQZuIqI0w8BNRJRmGLiJiNIMAzcRUZqxpByQiCiW0nIvFq6uwX6fH33dLsyaNhIlRZ5kDyttMHATka1Ky714YPkW+ANBAIDX58cDy7cAAIO3QUyVEJGtFq6uaQvaCn8giIWra5I0ovTDwE1Ettrv85t6nKIxcBORrfq6XaYep2gM3ERkq1nTRsLldLR7zOV0YNa0kUkaUfrhzUkispVyA9JoVQkrUKIxcBOR7UqKPIaCLytQ1DFVQkQpS6sC5Z4lFZi8YC1Ky71JGllyccZNRCkrVqVJR2ffiUrBSCkhhDD9vI7gjJuIUpZepUm89d9KCsbr80Pi1JuA2Rn8uprDuOr3H+Lw8QbTY+gIBm4iSllqFSiR4qn/7ugioAN1fty9eCO+99dP4G8K4qsTTabH0BFMlRBRygqvQPFqBOh46r/jXQQUCIbw93/vwtNvb0dzSGLWtJH44YVD0Cnb3jkwAzcRJZVerlmpQImsMAHir//u63apvhHEehMo23UEc0orUX3wOC4u7Inzh56J59bvwq9X19hepsjATURJY6bcz2z9dyyzpo00/CZw5GQTFry5DUvL9qFvfi7+9J3xqG9sxn+/Upm0MkUGbiJKmli5ZrUAaLT+OxZlhu8PBOEQAkEp4VF5EwiFJF7auBfz36zGiYZm3HnREPz84uHI65SNyQvWmhp3ojFwE1HS2NVwSgnWXp8fAoBsfTwoJVxOB6YUFmDh6hrcu6QCfd0u3DxhAN6pPoyNu49iaEFndHJkYdF7O/HapgOYNW1k0htlMXATUdLEk2vWopUrj0zHyIjn+QNBvLBhT9vjXp8fT6yuQeccB246bwCWf7oPDc2htp89sHwL3HlOHK0PJGTc8WA5IBElzZTCAkQuXYnnhmOsumy1dEykyGAOAF1znHhve21b0Fb4A0FIiaQ2yuKMm4iSorTci2UbvVFB89wB+VF5Yr3KE61c+cylmxCUamFZ38Fj2otq6vwBPP3tc5LW/IqBm4gspRV0tWbC/95xBKXl3rYgaKTyRKvGOyhlu5y2GUraQ+3cEi1vFsnqVGg4cAshHADKAHillFdaNyQiyhSxgq7WjTwlKIaX/2nNphWxgrNe0HYIIBhxUHjaI7JsUOH1+XHvkgqU7T6Cx0rG6rxKYpmZcf8cwDYAXS0aCxFlmFjlflo3JoH2QV0rwAelxAPLtyDXmWVqRu10CASCEgVdclDnD6Ap2D6H7XY5MXfG6KhUjNbMe/GGPQBga/A2dHNSCNEPwBUA/mztcIgok8Qqm5s1bWTUjUlFeHVGrEoNfyCoWt0RSyAocd+lI+B0iKigDQCdc7KjVm6unz1Vc6wA8MKGPba2mDVaVfIbAL8AEH2VREQaYu0vWVLkwc0TB+hWlRhpNGXWzy4ejgN16jcftd5sYr2BKOkdu+gGbiHElQAOSyk36hx3hxCiTAhRVltbm7ABElH60ttf8rGSsXj62+fA43ZBAPC4XZh/zdioGe/8a8bCkaCe157WAGx20+JYnxAAe3epN5LjngxghhDicgC5ALoKIRZLKW8JP0hKuQjAIgAoLi6Or/6GiDJKSZEHZbuP4MWP9yIoJRxC4NrxHtUmUnrnAbRvFBoV/qZhpl+JIq+TAyeb1F/fzl3qdWfcUsoHpJT9pJSDANwAYG1k0CYiUqPUaiu11EEpsWyjN658sDLzNipLAFeM7a05m1fOF2u2H34dDyzfohm07d6lnnXcRGQZs02k9JQUeQwtqnG7nHj4qlG45tx+7R4vLfdi8oK17WrK18+eGtd1KNQaVFnNVOCWUr4L4F1LRkJEGceKZkw3TujfVoIX6Yav9cf90wvRrXOnqJ91ZMd4vb0vZy7dZGs9N3uVEJFlzN4ANOKxkrG4ZeKAdsHLnefEsrsnYcG1Z6sGbaBj25XpjTcoJRZv2IM5pVt0z5UIDNxEZBm9qpJ4NDYH0btrLpzZWejcyYE5V5yFsgcvwfiB3WM+ryOzf6MliS9+vFf3mERgjpuILJPIXWsA4IlV1Vj0/k40hyRynVmYNW0kbps82NBz420hq7bxgpZ4G1qZxcBNRJaKVe6n1/VPcfh4A+5+/lNs3HO07bGGQAiPr6qBO6+ToTeCeMr/IvPiysYLDYGg6jL7BJWa62LgJqKkMHKzMBiSWLxhN369ugbHG5ujzmGmQiWe2b9WXlyLy6bd3hm4ichSZtq6hgfiTXt9mFNaiS3eOlw4vAc++OxL1fNH5qhjzeLNzv7NVr/4A/Z0BWHgJiLLxNPW1evzY/TDq3CyKYgsAdw6aSDmzhiNCx5fp5ujjrfkT+t5WluUaeW6uXUZEaU9vbauWpQViiEJLC3bh1cr9huqUIm35E/reVpblN04oX9Sty5j4CYiy+i1dTVSYheePtFboh5vyZ/Wz+v8AdXXfKxkrOHl8lZgqoSILBOrBK+kyIOm5hAefW0rTjQ2QwhAq5pOCax6DaniLfnTG6faaxppjmUVzriJkkDpmTF49uuYvGCtrU347RQrvbFm6yE8885nONHYjOvG90PZg5e0tVyNZDR3HO+CHysWClmJM24im3WkZ0a6USvBu/2CwXht8wGs2XYII3p1wdI7J+G8wS2rHuOptdZ7PSMLfhK9UMhqQlqw0qe4uFiWlZUl/LxEmWDygrWqH8s9bpehTnXpqqk5hD9/uBO/feczCAjcc8lwfP+CwXA62n/wDy/Ly3c5IQTgqw+kfDDtKCHERillsZFjOeMmspkVHfNS3YadX+Gh0kp8dvgEpo3uhYevGq2ZFlFyx6fTJxOzGLiJbBbvDbR09OWJRvzqjW1Y/qkX/bq58NxtxZha2MvQcxPdyzuTMHAT2ayjedx0EAxJvPifPXhiVTX8gSB+MmUYfjxlGFydjG/6ezp+MjGKgZvIZul2I8ysSm8dHiytxKa9Ppw/9Ew8evUYDOvZxfR5TqdPJmYxcBMlQTJrgK1yrCGAJ1fX4PkNu9G9cw6eueEczBjXFyLOlnmnwyeTeDFwE2UIoy1SE01KiRWb9uOx17fhqxON+M7EgbjvspHIdzk7dN5M/2TSEQzcRBkgWRUYO2pP4OFXK7H+869wdr98PHfr1zC2X37Czp+Jn0wSgYGbyAZWz4bNVGAkYiwNgSCeXfc5/vTeTuQ4s/DLkjG46bwBcGTZtJPAaY6Bm8hidsyGjVZgJGIs66oP4+EVldh7xI9rijx44PKzUHBGTgdGT2YxcBNZLFH1yLFmykYrMDoylv0+P+atrMLqqkMY1rMLXvzhREwaeqbh8VPiMHATWSwR9ch6M2WjFRjxjCUQDOGv67/Ab9Z8hpCU+MX0kfjBBUPQyaZtuigaAzeRxRJRj6w1U565dBMA4xUYZsfyya4jmPNKJWoOHcfovl3x5fFGLFxVgxc27MGUwgKsq65lxUcSMHATWSwR9chaM+KglO1m3nqB0+hYvjrRiAVvVuOljfvgcbvwgwsG44WP97Sb8S/esKftePYRsRcDN5HFwmfDXp8fDiHabadlZPPaLI09DoH2OWq9ihG9mXkoJLGkbC8eX1WNEw3NuPsbQ/HTqcNw6VPvx9zdPHIcZC0GbiIbKMHMaEVHZE5bK2gr9vv8qnnwe5ZUYObSCgRbn94tz4lHrhqt2j62an8d5pRWonyPDxMGd8djJWMwvNcZbec3gn1E7MHATWQTMxUdasfG0tft0nxOMCzmH60PYNbL7fPiJxqb8dRb2/G3f3+Bbnmd8OT145AlgNv++knbrDzf5YTPH73beaSOrpYkYxi4iWxipqLD7Mx1SmEBXgjLOccSCEosXF2Dq8/pi7krqvD8ht0ISaBzJwdmXjYCjiwRNXN3OgScWQKBUOyZf5xtScgk3cAthMgF8D6AnNbjX5ZSPmL1wIgyjZmKDq1jtayrrjX1HK/Pj+nPfICag8fbHjvZFMSDr1RCLTQHghLd8pzI65SN/T6/6jFAy041ZD0jhZiNAKZKKccBOAfAdCHERGuHRZR5zGxIq3ZsLPt9ftNd87aHBW1FrPm0rz6A9bOn4osFV2juXsNUiT10A7dscaL1W2frn8RvVEl0Gsh1nvov53Y5Mf+asapVGCVFHsy/Ziw8bheMZB/6ul0oKfJg8tDuhsZx7gC36f/E4Z8MZk0bCadKXxKfP4BBGb5zfSowtPRJCOEQQlQAOAzgbSnlxyrH3CGEKBNClNXW1iZ6nERpTan4OBqWSmhsDsV8TkmRR3eGC7Sftb/ww0m4ZeIAOGIkm7vkOPDdSYMMvSGovYYyti652plWpWKGwdsahgK3lDIopTwHQD8A5wkhxqgcs0hKWSylLC4oKEj0OInSWqyKEiO0Uifd8qJn7Y+VjMWO+ZfjyevHqc6Ks7OyMHdFleEZt8ftUv1koJfPNnN9ZI6pqhIppU8I8S6A6QAqLRkRUQbqaL+SWAtnSsu9mLxgbbvH+3fPw58//AKBkER2lkBzWDWIkbI+oGWWrZXKAYzdQGVdtzV0Z9xCiAIhhLv1axeASwBUWz0wokzizlO/aaf1eCS1FZEAUPToW7hnSQW8rZUeXp8fM5duwrV//Dd89U34n1vORU8TLVeVFEv46k6tdIeRG6jcH9IaRmbcfQD8XQjhQEugXyqlfM3aYRFlFq2FjzoLIgEAc0q34IUNe9pSG16fH7Ne2gSIljK9SEEp0SUnG2vuuwidc7Jx9+JPDY3R5XTg2vEeLNvoVV3dCbSf8U8pLECuM0tzoRD3h7SObuCWUm4GUGTDWIgyVp1GekLrcUVpubdd0FboLYQ52diMzjkt/721UhrhddnKLF4rF3/vkop2Y4hsMgUATodA507ZqPMH2C3QYlw5SWSDeFu7LlxdE1ftbWTpnlpHwEeuGh0VWO9dUqF6PiNjCAQlOudko+KRy+IYMZnBTuhENjCz+CZcPDf3nFkiqnRPqQkHYuevO5qT5s1IezBwE9kgckGNVoldpHgCaZfc7KjzKrvkuJyOtk6DarXWZldsJmK8ZB4DN5FNwhfUKPnkwTqrDOMJpFr11UZqySPfYGIt5NEaL1mPOW4im5nZab2kyIPth45j0fs70RySENDPN2vNeo3WkofvpBM51ljcLidvRtqEM24imxldRXn4WAN+9mI5/vDuDvTr5sLfv3+e7rlj5c21AnqsxlBqKZ5bJg5QzdfPnTFad3yUGJxxE9lMa7Wh8ngwJPH8R7vw5Fvb0RgM4Z5LhuOui4Yi1+mIuVrRo1OCN2vaSMx6aVNUKeHJpmaUlns1n6e2l2XxwO66GxOTdRi4iWzm0Ng/0iEEKvb68OArW1C1/xguHN4Dj149BoN7dG47Rqu0z8iNzpIiD+atrGrX6Ao4tbGCmcBrZGNisg4DN5HNtPaPDEqJb/1hPXqekYNnbzoXl4/tDRFxc1Bvs189WjcuWcaXXhi4iWzmiZHu+N75g3HvpcNxRm7svHO8s914FwJRauHNSSKbqZX4ZQng/102Ag9fNSpm0LbitdlTJP0wcBPZ7LLRvXD+0DPbvne7nPj19ePwk6nDLX/teBcCUWphqoTIRm9VHcS8lVvh9fnxX8X9cP/0QpzZxXjb1UTgjcX0x8BNZIO9R+oxb2UV1mw7jJG9zsBLd03C1wYZ2x+SKBIDN5GFmppD+N8PduJ3az9DlhD478sL8b3Jg+F0MEtJ8WPgJrLIv3d8iYdKK7Gj9iSmj+6Nh68axeoNSggGbjotqW0Flqi8b+3xRvzqjW14pdyL/t1d+OttX8OUwp4JOTcRwMBNGSpWYDbT5MmMYEjinx/vxhOra9AQCOKnU4fhx1OGIbcDbVKJ1DBwU8bRC8xaTZ7mrayKO3Bv3ufDnNJKbN5Xh8nDzsSjV4/B0IIuHbsQIg0M3JRxYnXfKynyaC7vPlofQNGjb8FXb3zPxDp/AE++VYPnN+xGjy45eOaGczBjXN+opepEicTATRlHr+90rA57SgMmvfSJlBKvVuzHY69vw5GTjbh10iDcd9kIdLVw1SORgjVJlHG0KjeUx40u71brkQ0Af3pvBwofWoV7llSgzt+Eey8ZgbkzRjNok20446aMo9X6dEphASYvWIv9Pj+EADSa9LUTPnv3NwXxsxfL8fa2Q22PBYISf3h3B/p3z+NqRLINAzdlHLXWp4POdOGFDXtObftlIGgDgDuvZRb9zrZDeGRFFfYdjU6xhOfPiezAwE0ZKXLfxHuXVKjGaocQCEmJfJcTxxubEYzYHea4P4CrfvchtnjrMKyndpUI+1mTnZjjprRVWu7F5AVrdXdKX7i6RnOCrWxq0DknG50c0ZUgzRKo9Nbh/umFeONnF8Kjkz8nsgMDN6UlpVbb6/NDoqUKZNZLm1D06FtRgVxvNqw83x8Iaf787m8MRafsLPazppQgpJE7NCYVFxfLsrKyhJ+XSDF5wVrNkr5IWns8GqVswqvkzN15TkjZUsOt1HsD8W8nRgQAQoiNUspiI8cyx01pyUxOuSNBWwCYUljQrkrlaH0ALqcDT3/7HJQUeVBa7sWslzchEGx5Ha/Pj1kvbwLQsSX0RFp0UyVCiP5CiHVCiG1CiCohxM/tGBhRLHbllCWAddW1misxAWDeyqq2oK0IBCXmrayyZYx0+jGS424GMFNKeRaAiQB+LIQYZe2wKFMZvaGoRy3XbAWP26W7EvOoxs7pWo8TdZRu4JZSHpBSftr69XEA2wDw8x+ZpnZD8YHlW+IK3pF7JxppDZJlsn2IctNRbyUmkd1M5biFEIMAFAH42IrBUGbTa/6kRatFa2StdnieOZLL6Yh67ViyBNptoqu2ElO5Kel2OeHzR8+u3S4ugSdrGA7cQoguAJYBuEdKeUzl53cAuAMABgwYkLABUubQSzmUlnsxb2VVW4rB7XLiynF9sGyjV7d3duRqyXyXE0IAvvoA+uTnYownH29tPQSjwtfhqK3EDK8auXJcHyzesCfqHFeO62P49YjMMFQOKIRwAngNwGop5VN6x7MckNSYKeHT43G7sH72VN3j3q05jEdWVGH3V/UYP6AbNu45mvDX0Louo88nAsyVAxqpKhEA/gJgm5GgTaQlkTcU9coBD9T58aMXNuK2v34CR5bAP38wAct+dD5umTgARlPdRksO9T5JECWakVTJZADfAbBFCFHR+th/SynfsG5YlImU1MLMpZs6VFsNaN8YbA6G8Ld/78LTb29Hc0hi1rSR+MGFg5GT3fKGUTywO17ffKAtHeNyZmmumDR681GrvzdvXpJVdAO3lPJDwPAkhUiTcpOxo0Fba4l52a4jmFNaieqDxzG1sCfmzRiN/t3z2r1+5E3GBo2gLWC8b7dWG1kugyercOUk2UItaMYrvNoDAI6cbMLjb1ZjSdle9M3PxZ++Mx6XjeoVtX2YWlWL1luIhPFVj3o3L4kSjYGbbKEWNOPhcbvaAmIoJPHSxr1Y8GY1jjc0486vD8HPLh6Ozjnq/6zN5JzVugDG2jk+vDSRyGoM3GSLRNyoU9IPpeVe/OqNbTh8vBEAMKSgM/51xySM7H1GzOdr5aIF2s+81dIcejvHE9mJgZtsYTRoRlI2OlBmuA2BIOaUVqI5rNB6/1E/th04FhW4I2fIUwoL2tWEAy1B+trxHqyrrm03kwbQts1ZX7cLJxub41o8RGQFBm6yhdoNPL2g7XI62vLZUkq8WXkQ9y2tQMQmNWhoDkUF0DmlW9ptVeb1+bFso1c1SEcGXrXZtZZE1aUTmcHATbZQu4EXK+h5woLq7q9O4uFXq/De9lrN470+PyYvWNs2W263v2QrfyCIddW1uotizOTjHUaapBAlGAM32SbyBp7eisPG5iCeWfMZnn33c3RyZOHhK0fhzx/sxP66BtXzK3nnXGeW5kzeSK7drl7fRPHi1mWUNLG2Afvgs1pM/80HeHrNdlw2qhfemXkRvn/BYPxiemHM1Zf+QDBmO1Uji2LMLJzR2oOSyEqccVPSqKVP7vj6EKzZdgivbT6AwT064/nbz8OFwwtUn2M2v2x0UY1aPt7pEIAEAmEJdi6yoWThnpOUEpqDIfzjo9146u3taAqG8ONvDMOdFw1BbozZtVaqxe1yorE5FHUj9OaJA/BYyVhD41Gr2Qa4yIasY6bJFAM3Jd2ne45iziuV2HrgGL4+ogCPzhiNQT066z5PbTWmUokCMMhSeuFmwWko1qq8TOWrb8Ljq2rwr0/2oNcZufjDzefim2N6Ry1V16K31DzT//7o9MXAnQJOt1V5Ukos+9SL+W9sg88fwO2TB+OeS0egi8ZS9Vi41JxORwzcKSDeLb3S0fZDxzHnlUr8Z9cRnDvAjedLxmJU367JHhZRWmHgTgGZ1IhfK+VzsrEZv33nM/zlwy/QJTcbj187FteP748sszv4EhEDdyrIlEb8aimf2cs2Y/M+H1ZVHsT+ugZ8u7g/7v9mIbp37pTk0RKlLy7ASQGxFqKkE7WUT0NzCM+t34WuLidevmsSHr/ubAZtog7ijDsFdKQRfypVo8RK7az86QVwOjhPIEoEBu4UEU91RKpVo2ilfDxuF4M2UQLxf1Mai1WNYrfDxxvQOz836vF0TPkQpTrOuFNEPCmPVKhGCYYkXvh4NxaurkFjIIRpo3thy746HKhrSHrqhihTMXCngHhTHmarUcLfHNx5TkgJ1PkDcQfYzft8ePCVSmzx1uGCYT3w6NWjMaSgi6lzEJF5TJWkgHhTHmaqUZQ3B6/PDwngaH0APn8AEqfeKErLvYbGW+cP4KHSSlz97HocOtaA391YhOdvP49Bm8gmnHGnAKMpD7V0yvxrxhpKsejt6mJkpaaUEqUVXvz/17fhyMkm3DppEO67bAS65joNXikRJQIDdwrQSnlkCYHSci9KijwoLfdi1kub2vpBe31+zHppExZeP053Ky6g4zu/fH74OOaUVmLDziM4p78bf/veeRjjydc9JxElHlMlKUAt5QG0bIulpDDmrqhq18QfaGnqP3dFlaHXMLIK050XPXP2NwXxxKpqfPOZD7DtwHH86ltjsfzu8xm0iZKIM+4UoKQnZi7dFLWHoZLC8PnVt+PSejyS2q4ukSJbs6/ZegiPrKiC1+fHtef2wwOXF6JHl5yo5xmpiEmlhUJE6Y6BO0WUFHlw75IK1Z8lorzPyJZfda1vAvuO1mPuiq1Ys+0QhvfsgiV3TMSEIWeqPsdIRUyqLRQiSndMlaQQrXRGX7cLeU71X5XW42pKijxYP3uq5ga3ffJz8cd3d+DSp97H+s+/xOxvFuKNn1+oGbQBYxUxc1dUpcxCIaJMwMCdQmKV9+Vo7L2o9bjZ1wGAA8ca8PiqakhI/GLaSNx10VDdpep6FTGl5V7NdE46tq0lSgW6gVsI8ZwQ4rAQotKOAZ3OSoo8mH/NWHjcLgi09PiYf81YlBR54KvXyHFrPK73OteO9yCyE7aS424IhPDE6hpDdd1anxIkWjbznbdS++ZpurWtJUoVRnLcfwPwewD/sHYoBGg3m0p0z+612w4j1jbRRnfgiXXTUyuXHv5cIjJPd8YtpXwfwBEbxkIxJLJnd6W3DvvrGnSPM5LKCP+UYEa3PCdvTBLFKWFVJUKIOwDcAQADBgxI1GkpTE52VtvMtlueE49cNdpU8DvWEMBTb23HPz7ahSwBhGJNuWF8Nq+M4R6NqphILqcDj1w12tCxRBQtYTcnpZSLpJTFUsrigoKCRJ2WcKqcLvwmX0MgZPj5Ukq8WuHFxU++h79/tAu3TByI68b3i8pxR6pvajaU51bGp8Xtcqrm7YkoPqzjTgMd2QV+R+0JPPxqJdZ//hXO7pePv9xajJ21J/HA8i1ROW5nFhD+fnC0PqBakx25kEat3E/hcjowd4a5TwZEFBsDdxow23e7tNyLJ1ZVt+Wxc51Z+OXVo3HThIFwZAncvfhT1UAbkgKA+spNpV9K5EKa8P4paji7Jko83cAthHgRwDcA9BBC7APwiJTyL1YP7HSlNqM1U1FSWu7FL17ejKZg2NRZAmfkOuHIakmOaAX8yOX2CuV4tZl/rKDtcbsYtIksYKSq5EYpZR8ppVNK2Y9B2zqRPbOVpeFTCgsMVZTs9/kxe3lE0EbLTuvhqxS1bjpmaSS9leZTZhfMsNyPyBpcOZlCtHLZ66prNRfmAEAgGMKi93fgkqfe07xpGR50tUoLc7LV/zkoE3EzNeMs9yOyDnPcKSRWLltrYU7ZriN48JVK1Bw6josLe6LSW4dDxxujjgsPuuENp8JTMlpNrpTmU0Y6DCpY7kdkHQbuFBJrQ4XBs19v1w71yMkmzH9jG17auA8etwuLvjMel47qhVcr9kcFV7W0itobgVbnQCXoRwb8LCFU8+KcbRNZi4E7hWjNaJXg6PX5MXvZZny04yus3noQJxqacddFQ/Gzi4chr1PLr1JrNm0kkKq9fmTQDw/4kVUmyvGcbRNZS0iNSoKOKC4ulmVlZQk/7+kgvKpEa0YLAOcN7o7HSsZgRK8zLHv9WEFfOc7r88PROk4PN0ggipsQYqOUstjIsZxxp5jwGe3g2a9rHrfkjokQQm/tY8deX0vkTFt5c6lvak74eIgoGqtKUlif/FzVxz1ulyVB2yitHeOVlZZGlskTUfwYuFPUri9PoqsrevPeeDsCJlKsem7ubENkPQZQbj60AAAMyklEQVTuFNMQCOLpt7fjst+8j31H/fhWkQd983NTqkGTXj03d7YhshZz3Cnk/e21ePjVSuz6qh5XjeuLh644Cz27qqdLkkmvnps72xBZi4E7BRysa8AvX9+K1zcfwJAenbH49gm4YHiPZA9LkzLjn7uiKmo/yVRI5RBlOgbuJGoOhvD3j3bjqbdqEAhJ3HfpCNx50RDkZJvfANhuSvWJ0fJBIkocBu4k2bj7KOaUVmLbgWP4xsgCzJsxGgPP7JzsYZlmpHyQiBKLgdtmR0824fFV1fjXJ3vRu2su/njzuZg+pndSy/uIKL0wcNskFJJ4+dN9WPBmNer8AfzwwsH4+SUj0CWHvwIiModRwwbVB4/hodJKfLLrKIoHdsNj3xqDwt5dkz0sIkpTDNwWOtnYjGfe+Qx/+fALdM3NxhPXno3rxvdDltaOBUREBjBwW0BKidVVBzFv5VYcqGvADV/rj/unF6Jb507JHhoRZQAG7gTb81U9HllRiXU1tSjsfQZ+f1MRxg/snuxhEVEGYeBOkMbmIBa9txO/X/c5srME5lxxFm47fxCyHewqQESJxcCdAOs//xIPlVZi55cnccXYPnjoylHordHZj4ioozIucM8p3YIXP96LoJRwCIEbJ/THYyVjE/46peVeLHizGgePNQA4tUN6xV4fNuz8iotSiMgyKRe4O7KEek7pFizesKft+6CUbd8nMngv37gP9y/bjEDo1O40ypdenx8PLN8CAAzeRGSJlErAKjureH1+SJwKgkYb87/48V5Tj8ejYq8P9y9vH7QjsSc1EVkppWbcajurKEHQyOxVa39Grce1qM36p4zsiSdWV+Of/9kDI6djT2oiskpKBW6tYBf5eHhgdec5ISVQF9FeNJzDRB+QyP0UvT4/Zr28CTnZDtQ3NeN75w/Gm5UHcKCuIeZ5soRo+6TA7nlElEgpFbj7ul3wqgTv8Mb8kYH1aL12wFbcOKG/odcvLfdi5tJNUTP0QFACCGLlTy/A6L75OLtffsyNBICWWf6slzcBEm1pFea/iSgRUipwq+2sEtmYX2ujWjXhVSXKLN3r88MhBIJSwhM2A1beELTSKoGgxOi++QBOBV1lJp3vcuJYQwCRae+WgN+ePxDEPUsq8OArW+B0ZKHOH+BMnIhMSanAHRkQ1QKakdyxAPDFgivavo+cpSvBOXwGrPeG4InYjiu8D3VpuRf3LKkwcIWnnGwKAjiVjuFMnIiMMhS4hRDTATwDwAHgz1LKBVYNSK8xv1Y6JfKYcHNXVGkGZX8gqLoFV7hY23EpbwodpYyDgZuI9OiWAwohHACeBfBNAKMA3CiEGGX1wErLvZi8YC0Gz34dkxesbbvRN2vaSLic2lt7RQbZ0nJvzKAMIObPHULg2vEeLFxdEzUWIPZM3ZElTNVb+vwBw6WPRHT6MhJXzgPwuZRyp5SyCcC/AFxt5aBi1XOXFHkw/5qx8LhdEAC65Tnhdjkh0JLOmH/N2HazVqP11M4sgZzs6L8OZRFP+FjuXVKBOaUts+xYqZtODgGHw1wLV9Z/E5EeI6kSD4DwFSz7AEywZjgttOq571vakkc2s8+h0Xrq5pDE09ePa7uBGYsEsHjDHry26QBilXT7AyFDrx2O9d9EpMfIjFttyhgVr4QQdwghyoQQZbW1tR0alFbwCklg1subTKUTIvPdsY4rKfJg/eypUTciteilYOJhdLxEdPoyErj3AQgvhO4HYH/kQVLKRVLKYillcUFBQYcGFSt4BYIS81ZWGT6XXk4cOJUXn1O6BUMfeEN3xm2Ey+lAtzyn5s+dDgFnxE44sW6CEhEpjATuTwAMF0IMFkJ0AnADgBVWDkoveB2tN34Tb0phTxQP7Nb2fbc8J26e0L8tR67kxct2H8HiDXtML4+PFH7OR64arfqm0S3PiYXXjcPC68dFjYNVJUSkRzfHLaVsFkL8BMBqtJQDPielND7ljUNJkUe3RE+tf0n4Uvg++bm4+KxeeLPyII6cbMStkwZi5rSR6JqrPgueuXRTh8ftcbuwfvZU1bFq1aUzUBORWYbquKWUbwB4w+KxtDN3xmjMemmTZhc+tf4l4Yts9tc14PkNuzGgex7+etsFGNsvP+brdXSmrZXmMHMjlYjIiJRq6xqupMiDb5+n3WMkMg+uVU/dHAzpBm3AXCOqSAKnuhiyDpuIrJZSS97DlZZ7sWyjehBUm91q3VDU6+KnpFc6MuNWnsml60Rkh5SdcWvNoB1CtLuJ5/X5cefzZZrniVWhEr7QR+11Jg/t3u7moRHcRIGIrJayM26tWm5lZhwIhvDch1/gN2s+g4TElWf3wZqth9DQfGrRi155ndabg9ZNxskL1hoqFeQiGiKyUsrOuGPNlH/x8mZc+Pg6zH+zGpOH9cDb916E3990LhZce7ap8jqjGzcojNSE642diKijUnbGrdabW9EUDKH2eCP+97vFuHRUr7bHzVZwGNm4IVJOdlbbmHKys9DYHL2sfUphxxYgERHFkrIzbqWZlJaglO2CdjzUZtBa6RUlHx5eW96kErQBYF11x5b8ExHFkrKBG2gJ3gVdclR/ZvRmod75wzsNxkqvqOXDtepQmOMmIiulbKrkeEMAT729HV+ebIz6WSJ7ehhNr5gJxsxxE5GVUjJwl+06gh+98ClqTzTi5gkDMKpPVzy7bkdSd0rXyocLtJ95s1EUEVktJQN3v255GNyjMxZ9txjn9HcDAG6aMDCpY9LayPja8R6sq65N6psKEZ1eUjJw987PxZI7JyV7GO0Y2ciYiMgOKRm4zQjvCGh1MGXDKCJKBWkduCM7ArJXCBGdDlK6HFCP1t6U7BVCRJksrQO32SXrRESZIK0Dt1a9NOuoiSiTpXXgNrNknYgoU6T1zUmW6BHR6SitAzfAEj0iOv2kdaqEiOh0xMBNRJRmGLiJiNIMAzcRUZph4CYiSjMM3EREaUZIqbUBVwdOKkQtgN0dOEUPAF8maDjpgNeb+U63a+b1mjdQSmlop3FLAndHCSHKpJTFyR6HXXi9me90u2Zer7WYKiEiSjMM3EREaSZVA/eiZA/AZrzezHe6XTOv10IpmeMmIiJtqTrjJiIiDUkN3EKI6UKIGiHE50KI2So/zxFCLGn9+cdCiEH2jzJxDFzvfUKIrUKIzUKId4QQA5MxzkTRu96w464TQkghRFpXIRi5XiHEf7X+jquEEP+0e4yJZODf8wAhxDohRHnrv+nLkzHORBFCPCeEOCyEqNT4uRBC/Lb172OzEOJcywYjpUzKHwAOADsADAHQCcAmAKMijvkRgP9p/foGAEuSNV6brncKgLzWr+/O9OttPe4MAO8D2ACgONnjtvj3OxxAOYBurd/3TPa4Lb7eRQDubv16FIBdyR53B6/56wDOBVCp8fPLAbwJQACYCOBjq8aSzBn3eQA+l1LulFI2AfgXgKsjjrkawN9bv34ZwMVCCGHjGBNJ93qllOuklPWt324A0M/mMSaSkd8vAPwSwBMAGuwcnAWMXO8PATwrpTwKAFLKwzaPMZGMXK8E0LX163wA+20cX8JJKd8HcCTGIVcD+IdssQGAWwjRx4qxJDNwewDsDft+X+tjqsdIKZsB1AE405bRJZ6R6w13O1revdOV7vUKIYoA9JdSvmbnwCxi5Pc7AsAIIcR6IcQGIcR020aXeEaudy6AW4QQ+wC8AeCn9gwtacz+H49bMnfAUZs5R5a4GDkmXRi+FiHELQCKAVxk6YisFfN6hRBZAJ4GcJtdA7KYkd9vNlrSJd9Ay6epD4QQY6SUPovHZgUj13sjgL9JKZ8UQkwC8Hzr9YasH15S2Bavkjnj3gegf9j3/RD9UartGCFENlo+bsX6qJLKjFwvhBCXAHgQwAwpZaNNY7OC3vWeAWAMgHeFELvQkhNckcY3KI3+e35VShmQUn4BoAYtgTwdGbne2wEsBQAp5UcActHS0yNTGfo/ngjJDNyfABguhBgshOiElpuPKyKOWQHg1tavrwOwVrbeBUhDutfbmjr4E1qCdjrnPwGd65VS1kkpe0gpB0kpB6Elpz9DSlmWnOF2mJF/z6VouQENIUQPtKROdto6ysQxcr17AFwMAEKIs9ASuGttHaW9VgD4bmt1yUQAdVLKA5a8UpLv0l4OYDta7k4/2PrYo2j5Dwy0/KJfAvA5gP8AGJLM8dpwvWsAHAJQ0fpnRbLHbOX1Rhz7LtK4qsTg71cAeArAVgBbANyQ7DFbfL2jAKxHS8VJBYDLkj3mDl7viwAOAAigZXZ9O4C7ANwV9vt9tvXvY4uV/565cpKIKM1w5SQRUZph4CYiSjMM3EREaYaBm4gozTBwExGlGQZuIqI0w8BNRJRmGLiJiNLM/wGiCrDbs9QxRgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from numpy import *\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "#读取数据\n",
    "def loadDataSet(filename):\n",
    "    numFeat = len(open(filename).readline().split('\\t'))-1\n",
    "    dataMat = []; labelMat = []\n",
    "    fr = open(filename)\n",
    "    for line in fr.readlines():\n",
    "        lineArr = []\n",
    "        curLine = line.strip().split('\\t')\n",
    "        for i in range(numFeat):\n",
    "            lineArr.append(float(curLine[i]))\n",
    "        dataMat.append(lineArr)\n",
    "        labelMat.append(float(curLine[-1]))\n",
    "    return dataMat, labelMat\n",
    "\n",
    "#标准线性回归函数\n",
    "def standRegres(xArr, yArr):\n",
    "    xMat = mat(xArr)\n",
    "    yMat = mat(yArr).T\n",
    "    xTx = xMat.T * xMat\n",
    "    #判断行列式为零，则无法求逆\n",
    "    if linalg.det(xTx) == 0:\n",
    "        print('the matrix is singular, cannot do inverse')\n",
    "        return\n",
    "    ws = (xTx).I * (xMat.T*yMat)\n",
    "    return ws\n",
    "\n",
    "#拟合数据\n",
    "xArr, yArr = loadDataSet('ex0.txt')\n",
    "ws = standRegres(xArr, yArr)\n",
    "xMat = mat(xArr)\n",
    "yMat = mat(yArr)\n",
    "yHat = xMat*ws\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111)\n",
    "ax.scatter(xMat[:,1].flatten().A[0], yMat.T[:,0].flatten().A[0])\n",
    "xCopy = xMat.copy()\n",
    "xCopy.sort(0)\n",
    "\n",
    "yHat = xCopy*ws\n",
    "ax.plot(xCopy[:,1], yHat)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "计算$\\hat{y}$ 和 $y$ 之间的相关系数，用来判断预测值和实际值的匹配程度。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-12-11T12:17:34.835824Z",
     "start_time": "2018-12-11T12:17:34.818812Z"
    },
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1.        , 0.97223133],\n",
       "       [0.97223133, 1.        ]])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "yHat = xMat*ws\n",
    "corrcoef(yHat.T, yMat)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. 局部加权线性回归(LWLR)\n",
    "标准线性回归依据最小均方误差（MSE），则可能会导致欠拟合。因此可以在估计中加入偏差，降低预测的均方误差。  \n",
    "__局部加权线性回归__:给待遇测点附近每个点赋予一定的权重，然后继续进行最小均方差回归。设置一个权重矩阵$W$,给每个数据点赋予一个权重。可以使用__核__来给附近的点赋予更大的权值，常用的是高斯核：  \n",
    "$$ w(i,i) = exp(\\frac{|x^(i) - x|}{-2k^2})$$  \n",
    "权重$\\mathbf{W}$只含有对角元素，点x与x(i)越近，上述权重越大，k为一个用户自定义的参数。  \n",
    "随后计算线性模型的权重$\\hat{w}$:  \n",
    "$$ \\hat{w} = (\\mathbf{X^T}\\mathbf{W}\\mathbf{X})^{-1}\\mathbf{X}^T\\mathbf{W}y $$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-12-11T14:58:04.598777Z",
     "start_time": "2018-12-11T14:57:48.535477Z"
    },
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAD8CAYAAABXe05zAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3Xdc1dX/wPHXYYMsBWUp4ELBzEoT0VKzgZlm2dKy8etbfuvbtKFlg2jYTttpe28rs9RMzQVarhyg4EZREJC9uef3x4dxRca9eFn6fj4e93HX+Xw+53OB9z2cz3mfo7TWCCGEaD/sWrsCQgghrCOBWwgh2hkJ3EII0c5I4BZCiHZGArcQQrQzEriFEKKdkcAthBDtjARuIYRoZyRwCyFEO+PQHDv19fXVoaGhzbFrIYQ4JW3YsCFDa93ZkrLNErhDQ0NZv359c+xaCCFOSUqp/ZaWla4SIYRoZyRwCyFaT3w8jB5t3AuLNUtXiRBCWCQ2FhYvNh4vWtS6dWlHpMUthGg9MTEQHW3cW8K8hX4at9alxS2EaD1RUda1tM1b6GC71np8vLHvmBijTm2cxS1upZS9UmqTUmpBc1ZICHEaa6wVbd5Ct7a1Xtcxqh5PnWp8CcTGNk+9bcyarpL7gMTmqogQ4jQXHw+XXmoE0KlTGy9f1Vqv3UJurDulqtUeG3t8Cz46GiZMaDAA67g4Si4dQ87KNaTlFrM/s4AdR3LZ9tI7bN6yl8JnnmviyVvHoq4SpVRX4DLgOeCBZq2REOLU1VCXRGws5OQ0vH1VoM3OBm/vuvczdSqsW2ccKzzceAxG2dhYCq+4ikwnT47+9x4yS0xkdoogM3osGSUmchYvI7fDYHK/3Erh+nKKSisoLjNRXFZBcVkFRSVlmAbcBb9nw+9La47ZZyL0mcivIzzpf/KfUqMs7eOeDUwDPOoroJSaAkwBCA4OPvmaCSFOPQ2NIomJMQIywKxZdW9f1S2SnW3sJyUF0tJg5kyYMoXisgpSXTuRGjKAQ15dOBR6NofOvJpDYf1J/W4vRyNup2ifC0TcDKtyjX0FXwiJRbhXlOLVpRee5cV4hoTQxcMFV0d7nB3tcHW0x8XRHtf0w7j+tRznM8Jxnf8TLtnHcCkvwamiDIezziI45i3bf2Z1UI0tFqyUGguM0Vr/Tyk1EnhIaz22oW0GDRqkJXNSCHGc+PiaLpBZs5p0EVBrTVZBKXv/Wseer35ib0YB+918ONQpgEM9IsjILzmuvB3g5+VCkLcrgeUFdNm0Dp+RQ/GN6I3vwT34fPohPndPwWd4FC7r/7b8AuXo0cYXh5cXBAVBYSEcOwYvvQRTplh9XgBKqQ1a60EWlbUgcD8P3AiUAy6AJzBPaz25vm0kcAshTlAV7KKj6x4FYtaNUnLuYPZmFJCcls/ejAL2ZhSwJ6OAvUfzyS0ur97EQZvolptOV39vAvv1JqijK0HertX3/l4uONrXcymvsfrUpaqOEybAvHk1Qd7b2+jm8fKq+a/BStYE7ka7SrTWjwKPVu54JEaLu96gLYQQdarqCsnONgJgVBTlFSb2ZRaSnJbHzvcWkuRyNknf7mLvr8eoMNU0KgO9XOjeuQOXnxVId193erz+At3//JWukQNwWLSw/mM21Kde1e1S16iU+rarr6snKMgI3EFBVn4oTSPjuIUQzc5k0hwMG8DO4LNJSskk6eN4dq4rY8/RAkorTACogEiCO2UTFtKZ6PBgwj5+m97LfqP7wAhcF9YahWx/C+Tth5gnGz5wQ33qDY0hr2+7+oL9Bx/UBPoW0GhXSVNIV4kQFmhnSR+W0FpzOKeYnWl5Riv6SD7J6Xkkp+VTVFZRXS7I1Y7ewT708fOgt58Hffw86NXFHVcn+5qdmX8+0LTPqqmfcSv8bGzax90UEriFsEBT+ljbCK01R/NLSDqST1JaHsnpeew8YgTovJKaPujOHs6VwdmdPn4ehPl70LuLOx4ujtYdsB1/VpayaR+3EMKGzFtyDfWx2uoYDbUWLSyXmV9CUprRck5KyzOCdXoe2YVl1WW83RwJ8/PgirODCPP3IKyLO2F+HnTs4GSbc2quz6qdksAtREuq3XfalNZjYwHX0hn3apXLKSwjqTo451UH64z80upNPFwcCPPz4NIzAgjzM4Jzbz93Ors7o5Sy/lwsZe2cJqc4CdxCtKSqFmNVarUlfai1A3VVwF2/Hn79te7REtnZsHOnMUztpZegf//qMdR5L71Gcmg4ybc8ws7A80judy5JM/8kLbdmDHQHJ3t6+3kwqm8Xwvw8qm9+ns6otWsh9mHjOL26GxvMnQszZlQnwojmJX3cQrQGa/pszcvGxBgBeNs2KCiAyEhYu7bObQqX/cUun27sDO1Hcvd+JJU5kdQ5mFTPLtXFXBzt6N3Fg952RfRZu4ywq0bT28OeoJefRVV9UdT+4qir7l5ekJsLnp6Np62LOkkftxBtlXkCB1jWZ2vevxsba8y94eVV/XZxWQW7j+aTnJZfPZojaeh9pAy4C62M5BMnUzm9ju5n8JEkekeGEHZuP/r4edC1oyt269bCuEmQmQlHNxk7Ne9qqd31Uld/c9eukJBg3J+Co2XaGgncQrSkqgmQsrPrbinXJSoKYmIoe+BB9ppcSBr7fyQNu5ikfekkdQ5h3+O/Y6oM0A4Kenbx4Mxeflxd3cXhTnDSFhyeedf4wvj6JQiLgYjKoBobawRtH5/jg3FVd07tL5nK+hw3VM/Dw2j9z5olq9q0AAncQrQxFSbN/swC4yJhVSt6XTJ7zptOub3xJ2uXA6G9AwhL/pex2/8mrDCDPvsTCe3dFcf4uBN32nmoEUSrujnMZ9czb0FXtZDNy1Y9r5oitSpoV+1n1y4j8EdH1wT1qv2JZiGBW4iWVNUijYlBa83Bdz5i54dfkzTxVpKCw0lKy2fX0XxKy03VmwR3ciOsuz8XrVtOWMFRwv53Mz0uHGpMivT910ah4jzITIFegQ0fv/bselUXOOtqGdcOwOYtafP91G6tywiQZicXJ4VoAaXlJpLT80hIzSXhcG71fZ7ZhEmBXi5GFmFlkkoffyOb0M2pnvZV7YuW1vQrx8fDuHE1LWVLAm1dfdfSn20zkjkpRCvKKSoj0Sw4b0/NZVd6HmUVxt+aq6M9fQM8iAjwJGL3v/T96E1633sbnnfcbt2BTjZoStBtUyRwC2GpkwheWmsOZRed0Io+eKyoukxnD2cjQAd6Vt+H+nTA3q4Zk1VEuyTDAYUAy4Jy7VEe9WxTWm5iV3q+WYDOISE1t3puaKWgu28HzurmzfWRwfQL9CI8wIMuHi4tcabiNNNo4FZKdQM+A/wBEzBXa/16c1dMiJNWe1iaJSuwxMaSu3wViR3eJ+EB/+pWdHJafvX0oy6OdvT192Ssvz0Rq5YQcf14+l4cVX9ftBA2ZskKOAFAgNZ6o1LKA9gAXKG1TqhvG+kqES2iduu4sefmw9t8fNDz55NaZCLhtbkkuPuRMPRiEo4WkVJeE4B9OjgREehJv0Cv6u6O7r6VXR2nwYx1ouXYegWcw8Dhysd5SqlEIAioN3AL0SJqt6gbSG4pqzCxe+pjJDh3I6HYgYSO3Uj48RDZjm7Q/zqUNtF9bzpnFmUwcUscEel76BceTOdff6yZPCk+Hm4y+yKoPVxO5usQLcSq/+2UUqHA2cC65qiMOE001Pdc33t1vV5PokeevROJe7NIePsnEuz6kPD1DpIWZBtdHeFX4FxRRt+0vVyaHE/ExHFEeNnTd+ZjdMjOqtlJZw947GGj87pK7S+K2uOVp00z5umYNk0Ct2hWFgdupZQ78CNwv9Y6t473pwBTAIKDg21WQXEKMs+6q1LV51zfxcLsbON1qA6WesgQjnz7EwmpuWxfmkzCpKdJiMrggLMXzImHriPo5F9EvwAP/i88uKar4+evcXh9OpSVQWmCsb9nyoy5NqD+ro/GMgJbeN1BcfqyaDigUsoRWAAs1lq/1lh56eMW9TK/QAg1wbhqlrt+/YwAGhEB27dX9yOXDYlij18oCddPISGngoTNyST4BHOstOb3N9TH7bi+6IhAT7p41DFPdFXftI9PzbSoQ4YYdfH0rGlN11d/a/9bEMICNu3jVsZv/YdAoiVBW4gGVc1uV5Xtd+mlx08D6uFBvpMriV0jSIjbR8LlD5EQdi073f0oNQEbC3AyldM3v4jo0m1E3HotEZkH6PvOy7g/MQOizjH2Ex8P91fOwjdvXv1dLFWvmaWiN9hF09AkUZLqLVqIJaNKzgNWAVsxhgMCzNBa/17fNtLiFvUyC4R6yBCOLI8j8b3PSRhzLQmOXiTsSWdfQc08HR3dHGtGdVS2onvs2orDM083PD+0eavamrTu2mrvu6plXt882EI0ka1HlawGJM2ruZ3i/2aXlFeQnJZPomNXEu96jcStuSQuWWKsW9hjHOwoIsRHEdG9C1dVZRoGeuLv6XJiV4ff0OODcF19z+YrzVS1uJui9r7NW+ZCtBJJeW8rTqExwRl/rSHxrU9IHD+JRKdOJB7OZVd6PuUm43fNxdGOPn4ehAd4mt08rF/5W4hTiKS8t0dNncO4FVvq5RUm9mQUGBMqHc4l8XAeiYdzOZpXAr2ugO1F+HtmEh5grF1YFaSrE1iEEE0igbu9s2a1kcYyCxuQU1hG4hFjno7Ew7kkHsklKa1m3mhHe0XvLh4M792Z8NIsIn7+kr53/x+dRg6zxVkKIcxI4G4rzMc2V61MUjuY1pWZZ01L3YJMQ5NJsz+r0AjO1bc8DmXXzHjn6+5EeIAntwwNJTzA6PLo2dkdR3u7mmNNvuD4Y5/iffhCtCQJ3G1F7ZVJ4MQW9MMPGytpP/xwTeC2ZgharSCfZ+/EjqAIdtgFkfDun+zILWdnRhGF9k4A2Nspevh2YGBIRyYPCSEi0LP+Ge8aC8xNWWtRCFEnCdxtRVUANg+AtZmvpG0Fk0lzIKuQRPdQEh98g8TkXHasWkbK+dPhfKOM1+4cwosyuG7nBsK7dCD8pafo7eeOi6O9ZQepas2vX2/8R1B77LQQwmYkcLc19bWg4+OPX0m7HnnFZez8M57Er+eTeP5oEk2u7DySR2FpBQB2lfNGD+jqzUQ/CP/2I8IL0vB/PhalAmHqx5APpCRA1zqCbn3LV2Vng5eXMWZ6xgzjHmrORYbRCWEzErjbmvq6HMwzDqOiMJk0KccKq/ugqy4YpmRV9kV3uwDPXTmE9+7AtYO6VfdFh/l51LSizYcgDh1ac6x164yujbq6NOqa47pq7cLISKN/vq6x05JVKITNSOBua+oIjPnPzGTnmIkk+JzBjvNHk/jOGnYeyaPArBUd6tuBM7t6c92gboTnpBL+0ZsEPPoAamgDXRVNGYJovo150PbxOX5xApkdT4hmIwk4bYjWmsPL15D43hckjLmWREcvEtbvYJ+zd3UZDxcHYzy0f00CS5ifB65OFvZFN8aa0R9VLXYvL+jbt/5VZYQQjZLFgtsB8zUMq1YETzySa6SAVwrxcSPCqYzwjasIH38R4RecS5C364kp4K2l9pSrp0DWpxCtRTIn25jswtLqRWYTD+eRcDiXXel5lFUYX5rODnb09ffg0jP8iahsRffxN08Bv6T1Kt8QS0bCCCFsTgK3DVUPu6tOATeCdWpOcXWZzh7ORAR4MiKsc+W80R6E+nTAwTx5pb2RC49CtCiLArdSajTwOmAPfKC1fqFZa3UyWmDdP5NJk7p8DcnvfU7SuOtIyqsgeWcKuzz8KDSuF2Jvp+jZuQPndu9U3YoOD/Cks4dzs9RJCHH6sGQ+bnsgCbgYOAj8A0xqllXerU2Lrqu8r2/NKIeMDOvrYLbfgseeZG+PfuzLLGBfRgF7MgrYnZ7PrvT86hEdAF3K8gk7lExvh1L6mPKJmHwlYRcPNYbdSaq3EMICtu7jHgzs0lrvqdz5N8B4mmOVd0smTDIPhHWVnzmzpsVde5taq6EUlJSTml3EoewiUjdu5/Divzh07nkcTNzLvp7Xk/7bMWB19W78PV3o0bkD1wzqRu+CdMJ++IzeFw3F+9MPa461bh1kJcCYRcefU+31Favek4AuhLCSJYE7CEgxe34QiGyW2phPfj96dN1BzTxY1xqHrLWm+Jb/kHXtjWTll5KxM52szQlkvvchmTqAzEV7yPQYwpGPNpH6+zFyyo7/b8O+y0D8d2cQEBzC8M3xdE/4k+7jLiI0ojuhT8/AzVRWOeStH9APJl1g1LNqRRQw7s0v0k2YYKSBp6XBvn0151A1EkPm7hBCWMmSrpJrgGit9W2Vz28EBmut76lVznyV94H79++3ujKZ+SWs33+M4meeoyRhJ8XunhT7dKbYzZ3iS0ZT3CWA4kOpFG78l/yI/uS6epBfXE5eSZlxX1xePVl/bU7lpfg429Ep7SD+2ekEerkQeOO1BHq7EOTtSuC7s+ky+yUc7r8PXnzRyADMyTHGKA8ZUvNlUbXALBzfiq9vSJz5WOeq/S1cWDPpkiyBJYTAxuO4lVJRwFNa6+jK548CaK2fr2+bpvZxx+3O4Pr319X5noOuwMXFGRdHO1yd7PFwdsTdxQFPFwfcnR3wyM7EY90a3EeeT6fwXvi4O9OpgxO+ydvp9NqLuD/+KGro0Pr7nGv3jZuvNv7BB0ag3boVCgtrUrvNVwqHhhebrb1orfR9CyHMWBO40Vo3eMPoTtkDdAecgH+Bfg1tM3DgQN0UecVleuvBbJ2clqdThl+ij7p56TwPb13m3VHrOXMa3jg6Wmsw7i0VF2eUj4sz9u/jU3Mc8/eqREYax4iMNF738bH+mEIIUQdgvW4kHlfdLCsEYzBGluwGHmusfFMD93GqAmdVsGwsONYVaBtjTbCPizPqUhW0tT4x2AshRBPZPHBbe7NJ4K7SlIDcHPuuK8g3pZUvhBB1sCZwt/3Mycay8k6mr/gkVo+p9zUhhGhm7X+SKfM5pSXtWgjRTp1ek0xJq1cIcZpp/4FbJjgSQpxmmqWrRCl1FLA+Awd8gSZOMNJuyTmfHuScTw8nc84hWuvOlhRslsDdVEqp9Zb28Zwq5JxPD3LOp4eWOud2PAm0EEKcniRwCyFEO9PWAvfc1q5AK5BzPj3IOZ8eWuSc21QftxBCiMa1tRa3EEKIRkjgFkKIdqZVArdSarRSaqdSapdS6pE63ndWSn1b+f46pVRoy9fStiw45weUUglKqS1KqaVKqZDWqKctNXbOZuWuVkpppVS7HzpmyTkrpa6t/FlvV0p91dJ1tDULfreDlVLLlVKbKn+/x7RGPW1FKfWRUipdKbWtnveVUuqNys9ji1LqHJtXwtLZqGx1w1gpfjfQg5r5vSNqlfkf8F7l44nAty1dz1Y45wsAt8rHd54O51xZzgNYCawFBrV2vVvg59wb2AR0rHzepbXr3QLnPBe4s/JxBLCvtet9kuc8HDgH2FbP+2OAhYAChgDrbF2H1mhxVy8+rLUuBaoWHzY3Hvi08vEPwIVKKdWCdbS1Rs9Za71ca11Y+XQt0LWF62hrlvycAZ4BXgKKW7JyzcSSc74deFtrfQxAa53ewnW0NUvOWQOelY+9gNQWrJ/Naa1XAlkNFBkPfKYNawFvpVSALevQGoG7rsWHg+oro7UuB3IAnxapXfOw5JzN/QfjG7s9a/SclVJnA9201gtasmLNyJKfcxgQppRao5Raq5Qa3WK1ax6WnPNTwGSl1EHgd+AeTm3W/r1brTUmmaqr5Vx7TKIlZdoTi89HKTUZGASMaNYaNb8Gz1kpZQfMAm5pqQq1AEt+zg4Y3SUjMf6rWqWUOkNrnd3MdWsulpzzJOATrfWrlWvYfl55zqbmr16raPb41SzjuH19fXVoaKjN9yuEEKeqDRs2ZGgLJ5lqlhZ3aGgoLbaQghBCnAKUUhbPqCrjuIUQop2RwC2EaD3x8cbyg/HxrV2TdkUCtxCi9cTGGmvGxsZaVl4CPSCBWwjRmmJijIW+LV0z1jzQ2zKIt7MvBIsDt1LKvjJl9VQZcyuEaG1Va8ZGRRnPGwug5oHe2tZ6Q052Xy0c+K1pcd8HJDZXRYQQotEAah7orW2tVzEPslWPJ0ww9jVhQsMBuL4AbcsvEQtYNBxQKdUVuAx4DnigWWskhDh1xccbwS0mhl09z2DWkmSO5pfQ2d2ZLgVZTMg10d/LywigjWxfHcRrmzsXZsyAmTOhf//jypeWm0h54XX27jrG3nfms8fZm1znAfiuy6DL9LcYPftxei5ebOynrn1XBWioafXHxNR8eVj7JdJEFiXgKKV+AJ7HmBDoIa312IbKDxo0SMs4biHECUaPJmtlPLMnPcKXfgNwc7QnPMCTjPwSUtOzKVH2jNyzgSvss7j4s9dwczq+bXl43FUsSNP8OuhStnv4U6HssLdTuDnZ08HJATcne9x2bMetMI9OJfn4uDlSeCyX1NAwDncJ5lBBORWqpqOhk5PCKzeLTI9O5JYZsXB41i4uv3AAgYPOwMneDsfKm5ODwmnLv/DGG5TefS/F776H/Zo1uLg44dI9GJenn8Lz/KHY2zVtWiWl1AZt4ULDjQZupdRYYIzW+n9KqZHUE7iVUlOAKQDBwcED9++3eCy5EOI0oOPi+OzVr3il+wUUOrpwfWQI91/UGx93ZwByV67h/Xd+5Qf/Mzns4oWbkz1De/ri7GiH1pr03BI27D+GBs7MTWVowhqcOnpTUVRE4fBRFIaFU1BaTuGuveQn7SbLxYNML1/cdDmB/p0I2L2d0ISNdO/qQ+iLT9HdtwPebk7V9Ts69kq+PubC15HjOezk0aRz/PmuYZzVzbtJ29o6cD8P3AiUAy4Ys3zN01pPrm8baXELIWr7aPJ0nu46nPNz9/NkzI309qs7OJpMmr/3ZfHzpkP8vc+YhM9OKTo42TOqrx+XnxVI9+QtNSNLcnPBxwcyMmp2UrtLpa7X6nluevJJdoSeQW5xGWUVJsoqTJSW68p7ExpwevlFnDZuwNSnD8V2DhTnF1FcWsbl11+M7523NenzsSZwWzsP7UhgQWPlBg4cqIUQosqS7Ud06PQFesp/X9cVa9bUXSguTuvoaOO+MVVlp03T2sdH6zlzrK9UdLTWYNxbquq4c+YcX1cvL2NfXl7W16MSsF634fm4hRCnke2pOdz7zSb6d/Vi9g2DsHv66bpHbdQemdHQELuqsv/+a7S0p0yp++AN7aOhUSmNjR6ZN+/4YYxBQcffNzdLI7w1N2lxC2EBa1qY7dSRnCId+dyfesjMP3VaTlHDrdzan4c1ZevTlFZ1Q9vVd1wb/CyxosUtgVuI1tLUoNJOFJSU6cveWKkjnliotx/KMV5sSndIXFzTA2NLb3cSrAnczTIft1ycFMICdV1AO0VUmDR3fLGBpYlpfHDzIEb19Tu5HY4ebXRRREfXPb76FGDNxUnp4xaiJZn3ndZO926OY9iiXBO8sDCRJQlpPDk24uSDNjQ9S/IUJYFbiJZki9To+HgYMsS4WXKRz9Z1qSvgV75miovjrWXJvL9qLzdHhXDLsO7W7bs+zfUl105J4BaiJVW1HBubE8Nc7UAZGwvr1hm3OoJuzownWH7p9Xxd5sPSMy/gwNsfUhEXd2Kwt7QVW9fxawf8qVM5smY9kz/dyCt/JDFuQCBPjI2w4AMRTdEaiwULcdqaXeDD9yOn4Zh4mHPt+nL387MImW+0IitMmvzicnKLy8gpKiO3uIy84nJy3/yW3ExncucspGR7Ic7eZ+J1SQC+GYfxvXQi9nsy2ZtRwI4jefy9N4vEI7noM6+vOWgKOO9Lp0+f6xi271/Omv053UL6ETRgEKXf/8LRvBKOzV9JxeefoybfiAJcP/0I9//+B/fIQbg/+zzufyzBHoxWb2WgL3rsSdIyCjiSW8x23wG8eet9lDi68MKE/lx3bjeUalrqt2icXJwUooV8+88Bpv+4laGHtuMd2IWluiPl9g74ebqQW1xOfkl5o/twMlVQamdf53supjLO6eLK4LO6M7jwCCFPP8qRvFJ23X4fyZt2sAUPNgb2pdy+ae01N3twdnZEA+UV+oT6DkjdyayUP+kx84lT9qJrc7Lm4qS0uIVoAev3ZfH4z9s4/+BWPv7qMRw6dSR9z0E+WrOPjPwSPF0c8XBxwNPVEc/qe0c8XR2M+22bcJ/2EPZ5uZg8PMi5+TYyFi3l6KhoylbH0ePoAYJW/IFd9CXw0CIgDC5fQxAwEKpHsOTv+409ew6TMjyaw3c/gJODHZ1nv0ynuBXYe3qgX3wRk4biz74gP+o88tdvJm/MOPL9g8jfm0LJ2r+xG3wudgp8V6/E/7KL8Xe1I/TVZwkqyYVZs46fQe8UHQHS6iwdN2jNTcZxC1Hj0LFCPfCZP/TIl5fr7HfmHp+ibW0iifm477g4Y1+gdUSEZanfdaVs11eH2uPMzZ9XPY6MrKmDeb1O8cSi5oAk4AjRNhSWlOvL3lip+z25SCen5Z5YIDKyJgA2JC7OKBMZeWJmoY9PzX4sTeYx37a+AFs7ANeVEFN13Ib2IywigVuINsBkMum7v9qoQx9ZoJcmHqm7kKWBuy4nk1lo3lo/mcxNaV3bjDWBWy5OCtFM3vlrFy8t2sm00X3438hedRdqzezJUzhzsz2SzEkhLNVM2YNLE9N4efFOLh8QyJ0jetZfsDUTSySppd2SwC1OXZYE5alTjREQU6davk0jdqXncd83m+kX6MmLV50p45mFzTUauJVS3ZRSy5VSiUqp7Uqp+1qiYkKctLrmd24oVdx8m3HjmjTXR05hGbd9uh4XR3vm3jgIV6e6x1wLcTIsaXGXAw9qrcOBIcBdSinJZRWtr3bgrP28dkp3bCx63TpyN21l141TWDN/JfMeeIH3r7qPzwMHsujHFewbdw0mRyfIzDwxnbyR1O/isgru/nojh7KLeG/yOQR6u7bAhyBOR1ZfnFRK/QK8pbVeUl8ZuTgpGmSri2K1pvqsGBLFwZ37ODL4fNJemk36lh2kLV9D2tlDSHNyJz09m7TcYgrtnRrcrVdRHqMObuGiSdEMHzsUDxfHE46nFy4k+Y/VbP50HluHX8aWClcSDx6jVCte9D3GdQ/VuySrEHWy6WLBtXYxQTz5AAAgAElEQVQcCqwEztBa59Z6T1Z5F5YxD7gxMccFcT1nLkUxT1EWE0vp5JsoqzBRXqEp3bCR8nfepfy/d1B+Rn+KyyrI2LiV1F8Ws3PIKHaaXEg+kEmJveNxh3IuK8HfVIRfnx508XTGz9MFv0/n4rd9M13ys/A7OwKfB++h9IYbOVxuT2JwBOvCI1kefBbHSjUOdopgHzd6+HYgtDQX5zWrODBwGPE5kJFfCoC7swP9Aj0586fPGZWwmqjcA5Cd3RqfrGjHmiVwK6XcgRXAc1rreQ2VlRa3qFd8fPWFwPxb/sPiTxYQ79OdI6F9SO0RzpHUTAodXazaZRcPZ/r4e9DHlE/Y8gUEXnclfkPOpsvOLXjOfAZVu2U/dy5Mm2asD/jBB8Z7c+fC3XdDWRlER1P+2+9sPJDNiqR0dqcXsDejgH2ZBZSbNH4ezpzbvRPDevkyMKQj3X06YGenoF8/SEiAiAjYvt2Wn5o4Ddh8rhKllCPwI/BlY0FbiAZVTkm64urbmbqrI1kjp+BbkE23Tm709ffggrUL6XxgN06dO+H4+GM4vvE6jlv+xaFHKI6uLjhcdy2OduD0xef43n4z/ucPwcvNrJV9+4U1j/2GgWPMid0y8+ZBTo5xkdL8tbIy8PGBmBgc7O0Y3L0Tg4/shFnG9qbI0ej338f+sRkwcyYMqrVA7Qcf1BxLiObUWIYOoIDPgNmWZvVI5qSoj2nNGj37psd1yPQFOvq5hfrviCHa1NA8F3Vl5tWeQ6OuMrVTss2zAxsq39h8HVXZhj4+tvlAhKiELVPegfMADWwBNlfexjS0jQRuUZ+3lyfrkOkL9NRvN+mi0vKmpUxbshq4+SRIJ5OSXftYc+ZYNpmTEFayaeBuyk0CdxOcBnM+fLVuvw6ZvkDf+/VGXVFhst2OrWlBC9FGWRO4JXOyrbDl+n+trVadtNZ8tz6Fx37ayoiwzrx89QDjYp6t1JW6Lenc4hQmCym0FRMmwPr1xr012uKk9WZ12vLBNzz/+w7i92QytKcP704+BycHaS8IcTLkL6itmDfPyNab18CgHbOW7NG8EhZtO8yL1zzMomv/h+nJJxs/RmOZhrYSE0PKuGu497IHufytNexMyyP28n58eutg3JykrSDEyZK/oraiaghZ1erfdWQVbnswhk+9I1j/5U72umcBoBTo7mPotbaM/7kcZNyAQBztj/8+zi8pZ82uDP76cCUJXaIJ+PRvQrK96T73U0bErSdg6lRYu/bk6l+ZDZkz4wnezO3IZ/3/D7v0cu6+oBf/HdEDj03r4bI7ZApRIWxA5uNuQ5LS8nj3ibkcySuhk6cbPmMuwqeDMy6OdmzYf4w/EtLwKCkgMnUH5/7vBgaFdiIiwJMliWm8s3wXO47k0bWjK/8d0ZNzQzuyOjmD5TvT+XtvFmUVGncH6H8shfSAEFKKoLTChJ2pglEHNnP9dSMY4eeI/dNPNym46tGj+elgKc9F30mWoxtXZybwwKShBIwaZhQYMgTWrYPIyJP/khDiFCSLBbdDX607wOM/b8WtS1/6OqWyw78rmf+mkl1YBoCvs+KepD+5bfW3eD37FJjN8Xz5gEDGnRnA0sR03lq+iyd+3lb9XpifO7cO687IPl0YGNKxun+5wqTZ++ca5r32Bd/1Oo8/1+UTVJrLdXkdGff8bILmRVrcF11YWs7j4x5kXkop5+Sk8Gnyl5yx9BdIiYZRbaTfXYhTiLS424AtB7O58p04zuvly+yeZXR8/pnqVm9ZhYmisgo8rxh33IRK9dFas3bBKg58/TPD/jOBrheeV/+BK9PPS5Udfz70PF/tK2b10XIAFJrOHi4EersS6O1CoJcrAd6uBKYdIPCLDwm49794njeEv+av5qUV+9nr2on79q3knu9fxX7wueDtfXzLXVZbEaJBzTbJlKUkcFuuvMLE5W+tISO/hCUPjMDrynoCtDWBr9asedaU29cjgnUeXUn1C+bwlHtIzS4mNaeI1OwiistMde4m5FgqMxe/zTB/lxMDthDCItJV0o58vGYfCYdzefeGc/Byday5SFl1bx6wLR3uV3sfVpQLfeR+QmfMgLuuhqsHVL+utSa7sIxDK9dx+OOvSB1/LZklJgZ+/yFDk/7G8dhRmPEyTJlS+yhCCBuTFncrSskq5JJZKxnWy5f3bxpY9xJXlraeW0NV3Xx8jKGMbbGOQrQT0uJuB7TWPP7zNuwUPD2+X/3rElraem4N5kMY581rm3UU4hQkLe5WMm/jQR747l9ixkXwf8O6t3Z1hBCtzJoWt2ROtoJVyUd57KdtDArpyE1Roa1dHSFEO2NR4FZKjVZK7VRK7VJKPdJstbFFCvbcueDra9w3J/O6WlHvPxPS+M8n6wnxcePdyQOxt+VkS0KI00KjgVspZQ+8DVwKRACTmmuV98xnX+Rep/5kPfti/YUaC5gzZhgXymbMsO7g1gZi89n8qh5PnXridpX7MsXF8eOGg9zxxQbCAzz4ZsoQOns4W1dHIYTAsouTg4FdWus9AEqpb4DxQIKtK7PrjgdYvOoYO70u5Iu8kroDm/lseFD9uOL3haxMOspP971LevJ+irv3pHj2SorKKvCqKKH3nm2E9Qqk3/Jf6Z+XitfLLxw/1rie/dY7SqLyQpzpygkkf/crWy65kQR7T9KcPMj7fBP5m0zkF5eTf/Aw+eG3kv9LJlodY3DmHj4cPRCPLz4xvlxmzpQhdEIIqzR6cVIpdTUwWmt9W+XzG4FIrfXd9W1zMhcn43Zl8J+P19GpIJuYLT9zXvERyl94gayIs8gqLCVrwxaKvv0B+2uuxg4o+/4HVg8fz9IsTUZ+KT4dnOjZ2R0XJ3tcHe1wyTlG5rqNJHkGkO7hU30c/9J8ws7ojk8HJ2Ml8aMZlOxIJrd7LwrLNV57kuh8LI2gQf3pGuJHt8/fx6eskPz7HuBw1x6kbkli89rt/OPmT7ajGwCupcUEFWTi0SMYd9+OeLg44L5vNx3iV+NRmEdo+n7G7FiNy0WjjClcMzONoXQZGU36rIQQpw6bZk4qpa4BomsF7sFa63tqlZsCTAEIDg4euH///qbUHYBNE27mIZ8odvt0s6i8h7MDI3ztuGz1z1x47w04DRta86bZWOPsZ57n3/nLSXT3I2nEGHaaXMgrLsfBTuGQk41z6kE8gwNwDQ0hZ9kK0l08SfXsQpl93f+YhGalMrgkjcGlGZw97ExCF3yPfcyTx7fkq44fEQGHDtWsLL51q7S4hRDVrAnclqw5GQUsNnv+KPBoQ9uc9NJlcXG6+Iz++o9+w/WbwyfrOYOv1N9PvE8vS0zTmw8c00lHcnXi4Ry97VC23nYoW5eUVdS97mDlvupdwsr8vdqLwEZEaA26PKKfTl26Wv99ydV64Rkj9MrQs3TyBWN0/so1xnqGkZENL48lS2gJISyAjRcLdgD2AN0BJ+BfoF9D29hkzcmqQBwRYdnirE0JkObBvvYisHWtNu7lVbMArdY1K4hXPRdCiCayJnA3enFSa12ulLobWAzYAx9prbdb+V+A9aqy8LKzISHByMxrqEuhao3BphyjalIk8/3X3l9sLOTkGH3Ss2ZZdxwhhLChtp852VamA62rHm2lbkKIdu/0mtZVgqcQ4hRwek0y1RZXORdCiGbU/gN3W549TwghmkH7D9xNuSgphBDtWLP0cSuljgJNycDxBU63NEI559ODnPPp4WTOOURr3dmSgs0SuJtKKbXe0s75U4Wc8+lBzvn00FLnLPNxCyFEOyOBWwgh2pm2FribefWDNknO+fQg53x6aJFzblN93EIIIRrX1lrcQgghGtEqgbuxNSyVUs5KqW8r31+nlApt+VralgXn/IBSKkEptUUptVQpFdIa9bQlS9cqVUpdrZTSSql2PwLBknNWSl1b+bPerpT6qqXraGsW/G4HK6WWK6U2Vf5+j2mNetqKUuojpVS6UmpbPe8rpdQblZ/HFqXUOTavhKXTCNrqhjHD4G6gBzXTxEbUKvM/4L3KxxOBb1u6nq1wzhcAbpWP7zwdzrmynAewElgLDGrterfAz7k3sAnoWPm8S2vXuwXOeS5wZ+XjCGBfa9f7JM95OHAOsK2e98cACwEFDAHW2boOrdHirl7DUmtdClStYWluPPBp5eMfgAuVUu15OfRGz1lrvVxrXVj5dC3QtYXraGuW/JwBngFeAopbsnLNxJJzvh14W2t9DEBrnd7CdbQ1S85ZA56Vj72A1Basn81prVcCWQ0UGQ98pg1rAW+lVIAt69AagTsISDF7frDytTrLaK3LgRzAh/bLknM29x+Mb+z2rNFzVkqdDXTTWi9oyYo1I0t+zmFAmFJqjVJqrVJqdIvVrnlYcs5PAZOVUgeB34F7OLVZ+/dutdaYq6SulnPtoS2WlGlPLD4fpdRkYBAwollr1PwaPGellB0wC7ilpSrUAiz5OTtgdJeMxPivapVS6gytdXYz1625WHLOk4BPtNavKqWigM8rz9nU/NVrFc0ev1qjxX0QMF8FuCsn/utUXUYp5YDx71VD/5q0dZacM0qpi4DHgMu11iUtVLfm0tg5ewBnAH8ppfZh9AXOb+cXKC393f5Fa12mtd4L7MQI5O2VJef8H+A7AK11POCCMafHqcqiv/eT0SzjuH19fXVoaKjN9yuEEKeqDRs2ZGgLJ5lqlq6S0NBQWmwFHCGEOAUopSyeUVUScIQQrWbPktVkXnaFsQShsJgEbiFEqzhWUMr4P9KZ7nKmsQShsJgEbiFEq3hv5W7y7J1Z3utc0qc/YdlG8fEwevRp30KXwC2EaHFpucV8GrePwaGdqFB2/GTnb9mGVYuDx8baNoi3sy8EiwO3Usq+cq6BUyVZQgjRSt5atovyCs0r1wzgnGBvfthw0EgXbyyAxsRAdLRxbx7ET5IpNpbiP5c1fV8tHPitaXHfByQ2V0WEEKeHlKxCvv77ANed241gHzeuGdSN5PR8/j2Yw2dvz2Pnhh31BtCt3SKYee8sMvufc3wQN9dAENVa8+dPK1h43V0QH0/Rqjg+mTyNCwbdydD7v2LDmIkNB+D69m3DLxFLWDSOWynVFWPukOeAB7TWYxsqP2jQIC3DAYUQJ4iP58GP17DAN5wV0y/E38uF3OIyBj/3J4HOsCffRGjOERYNqMDljiknbH7BMwvZW2AiyNWO2y7qy4Szu+Ll5ggYXwg/bTrEknl/ce6Ov5lWlIDLq68YwTQmhn8C+vD874lsPGAkqT6bspwVniEs8epB36IMUjv6U1hUwpNL5nJTl3JYtOjE+o8ebQRo81Z/5RdH6dPP4PTkExAV1aSPRim1QVu4XqWl47hnA9Mwst2EEKJJkl94g5/6TOK2Ixvx97oMAE8XR0b38+fnzak4lZexz8ufN9f+w8N3HL9tRn4JewtMXJa4kg09ziL21wTenL+Z83ZvILXfOawvcQagH4qPzr2CjWlnMGf6E7BpE899Es/8jln4leYzM7CMpbuzebzbBQA8nBrPXXdeRk6ZiamfJ/DkJXfS73wPBpodW2vN4u1HeOm8+xgeNJyoy4bx2jdJ5PScRJ83l3FN3m5mDrmHOd0i6N8Cn2OjgVspNRZI11pvUEqNbKDcFGAKQHBwsM0qKIQ4NSz/eQX39r4aN1M5d9w06rj3ro8M4efNqUzft4zEDn7MCTyXsYdzCQ/wrC6zubKlfItDOs+sn8uKLPhs4Djm94qi15FDPDxxFFecHURQ4mYW3f0UD4y4ncsG3U5RpBNljk7cm/o3d375Iq4XjmTCr79x51tLydx7kNv/OwaiovAaPZrZy1dxzv1fs8SuS3Xg3p9ZQMz87fy18ygAezoP4JO/8wnGjqGp2/gjdBArgs8iICeHbp1cW+SztGTu2ecxcu/3AUeAQuCLhrYZOHCgFkKIKum5xTr8oXk6ZPoC/daNM+osczSvWJtMJp2VX6LPefoPfflbq3V5han6/ZcWJeqej/6mC0vKtY6L0zo6WpdOm653d4/QpvfmHL+zuDi9Y/wkfdGzC/WtH/+t92XkV2+j4+KMMtHRugJlvFa5jY6O1pNeXqQvfu0vXVRarmct2al7P/a7jnhiof5g1R6dmV+iV+xM1yuvuk0XODprHRmpDw6/WM+67E69MzRC6zm16mEFYL22dE5wSwsa+2UksKCxchK4hRDmHvnxX93rkQV6x/hJ2rRmTd2FzALrz5sO6pDpC/RHq/dUvz1xTrwe+8aqE8o2WT37eH/lbh0yfYEe9sJSHTJ9gb7ryw36cHbR8dvMmXP8tl5eRjj18mpydawJ3DKOWwjRrLan5vDNPyncPKw7fabfg3r66bpHbZiNzLh8QCAjujjy4i//8tpHS8kqKOXfg9mcE+x9QtkGNTRMLyrKuABZ62LiReF+2CtwPJzK51EevHX9Ofh7uRx/3Hnzjt82KOj4++ZmaYS35iYtbiEsYItWYxtnMpn0dXPi9Fmxi3V2YalxvlDTPWGu1ueRdtmV+o4rHtUh0xfoPo//rkOmL9A/bTxYZ9l6NXS8Buwdd60utnc4cbv6jmuDnyXN1VVi6U0CtxAWaGJQaU8WbTusQ6Yv0J/F7TVesCbAVZZNXLRS3/n6Yn32g9/rw0tXW1eBpgbUVvhStSZwN8t83DKOWwgLxMfXjANu4tjfFt2vlUrKK7hk1kqc7O1YeN/5ONifRM+s+fjpusZXnwKsGcctfdxCtCTzPtd6+lhP1qpXP+JQ3MaT6/+1gc/i9rM/s5AnxkacXNCG+rMkT1MSuIVoSZUXt3KefYGCknKLNikpr6C4rKLmhfh4GDLEuNUKuhv2H+PGXldw+e1vs/nexyyqi9Vp2nUF/FqvZeaX8MbSZEb17cLwMIsWdWlYM33JtVetsViwEKetDffM4MeuI/mh8xnwxG/c1Nudx2+/sM6yh3OK+Cx+P1+v2YN/7lF+vbY3GSWahz9eTX/nvtwd/y0dpk6FtWsBMJk0sb9up4uHM84VdkxcnkHU58/i2a8PnnYmui6ezw3HEnGf8zZERVH+xJMcc3DFd8bDda5uC1BeYWLdb6sZ/M4LOMZUpnNXBXyo7rb45JWvSLDrw/89/jzhS+fz2pIkisoqmDEm3NYfoaCZ1pyUPm4hTrTxwDEmvBMHwNWHN5NRWM5fPQex9alL8HBxPK7cx2v28fvWw2itGZR3iL/dg4g5uIJF3j1Z594VgMDcdOZt+wr/lUsA+O6fFKb9uIXZ153FsCfuJtY+jP3eAeS6eZLr7MYxFw8i0nbzXdoS3H//lZm/JzJ35R58nBRhaXsotXPgWFAIeWWaC1O3Me2GYbyc5cnXfx9gzI7VvF78L46Lfidu/kqW/LKKqTePxHP4MFKyCrnw+SWU2hvtwGG+DsQfLeWmnq48NeWiFv6U2y9r+rglcAvRArIKShn7xiryMrP57rOHCO/sxrJ+53Nrj3F8f0cUZ3XzZuG2I3y0ei+bU7LxcHFg4rnduCkqlK47NnPDl1uI8zAWDn/pwJ/0GDKA6w92JMzFxHNbfyLg4uGM2e1FiJ8nPzx8CWrtWrjtNjh0CF56CYBlb33F7aMf4NycFN6ZdBaXrC4kyNuVsH/jSS6xo0NpMd5eHbBDs9CzBw5KU2znSKSPA+syyxkd4MilI/px3zebAQju5Mbb15/DJ3H7WPDPXn755H7+vOBqPukWiS4tZeme7/D+7ZdW+8zbGwncQrQhJpPmlk/+Ye3uTH5cN4f+S3+ByEiOLP6LIc8v5fzeviSn5XMkt5juvh34v2GhXHVOVzo41/Rk7vpjNbE/bOTcnBTuue9K1NChLElI45EPV5Lp6AaAU3kp8zZ/yhl//lxvXX6ZeC8PdRuFuy7nmIMrc24cSHTuXpg61SgwaxYASS++ycyzryIoJZlnbh3OJ6YAnl6QAMCgTg78r7cLj606TKaLB2UmuP1AHDP2LIVZsyiu0JQ8/yJejz8ifdJWaI7ZAYUQTfTmsl2sTDrKzCv703/sdHAohpgY/Dyd6eysWJWcwfmdHXj+lnMZEdYZO7sTe5x7vfYsn1f1K+dshUWLuDhnD4PWvcOfvr1JsevAqLW/c0ZwxwbrMv6+SQS88j7/7XMFZ+ancXHOHhg6tLqfvErYz1F8UjUE72A0ty5aRId33mTTsXKeVHtwW1bObyvjefj6p9ji1oU7f3oTRgyFqChcABdpaTcvSwd8W3OTBBwhDCuT0nXoIwv01G82aZPJdML7O0aN1Uk+3bSOjGx4R3FxRpnIyOMmSapO4LEyYSR7zDid6+SqtY9P/dvU3qf586rHkZG6TNk1vB9hESRzUojWl5pdqM9++g998Wt/6YKSsroLRUYaf4aNBe66nEx2X1ycEWxPNnPzNEjbbynWBG4Zxy1EMygtN3HXlxspKavg3ckDcXOqp1dy1iwjsaSyb9kqJzO2OSoKfv315JNaZHx1q5DALU5dlmQGNlP24AsLd7DxQDYvXn0mPTu711+wNQOfBN12q9HArZTqppRarpRKVEptV0rd1xIVE+KkNZIZqLVmTezrZK6MrxlVYUUg3/jbKnIuu/yEsr9vPcxHa/Zyy9BQxp4ZeNKnIURtlrS4y4EHtdbhwBDgLqVURPNWSwgbMJvfQmvNtoWrSBkZjakyVfzdFbu54awbue2qJ6mozB3c9cIb3OF6DttuubvB4L3lYDYTVuVyYfdr+emNb4wLRsCeo/lM+2ELZwd7S9agaDaNBm6t9WGt9cbKx3lAItBCs4ULUb+KuDiSx0+qDrDH3v2Al0f/l11vfXRC2e83HGTsilzOH3IvE3tcwR93Pckri3YCsCmoLx92i+To8jXc3H8ii/oMY8aQyZhqtdSX/rSCP669E+Lj+WTNPjrYQ5AqYWrIJUycu5atC1dx57M/4qgrePv6c3BykJ5I0Tys+s1SSoUCZwPrmqMy4jRho37lV99fwsXhk1n1qhGon1uyi7fPupzoA748NX87E7/Ywm0dBpP+3Ms891si53R04OE9y/g7uD9TLrmfsIz9JFzWkUuydzOz58VcOP8QmSWaW9bPZ0tAGPNvfaT6WH/tTOf2tblMDRrFnhde59ctqVw9OISfXr+VmVf2Z8eRPMatyGWnqy+z//mCQO8WWjRWnJYszpxUSrkDK4DntNbz6njffJX3gfv377dlPcUp5JHbX2JvgYmzyWVsegJn5B8xMu4GDSbu7S/4Pn4vQWf1pe8l59G7izu/bE7ly/i9jMzazRPXnEPXC88z5sd4ZTmlJujrac/TkwZz7Zx4bkz4k4rIKL7K61B9vEBXO46Wahbedz69unjwy+W3saLCk8eWf4jP8CiOTn+C2z7fSJnWzFz9CWf6OHH5+KfIwpFlD43kQFYhV70Th7e9iZRCE6Ed7NhXYGLpgyOqLzxm5pcw67anCdifzF2522H79tb6eEU7ZfOUd6WUI7AAWKy1fq2x8pLyLuqz8bdVTFiVS7eiYxxx9qTMzp7wtD1kdexCmlPN6Avn8lJKHJwAsNMmRmUms8YzBO1gz13R/dh6KIdVyRk8FN2HZxYk4OHsgLuLA0sfHIGbkwOJh3NZkpDG5pRslu1I595RvXjgkj7GzufOhWnTjPUBP/jAGFUxdy7cfTeUlUF0NGvf+ZKJc9fyn/O6s3j7EUrKTfxy1zDu/3Yzf+/NYkRYZz69dfDxJzdkCKxbB5GRJ2QiCtEYm6a8K6UU8CGQaEnQFqIhn/7yD16ugSxK/hbTI4/yyWtfs9a9K/06KboNCaPby08zauXPuPcMZd+yeJLumU6vZQsI6xnAoc5dee7iKby6JAmAB/q6cuuwUH7bklo99K5qvHT4vu2EvxxLzvirWJCyg6vdzFLB582DnBwj0FYNhZs3zwjaPj4QE8OQHj5E9/Pjw9V7cTaV893ITgR6uzKlMJm/8eE/uYlArcA9a1bNyjNCNKfGMnSA8wANbAE2V97GNLSNZE6KupRXmPSAx3/TD9z+ct1p21o3nGZdaeVVt+nYUbfpwtGXaa21PrBklf5q0gPatGZNzcGq9ltXdmBd2X51vLb3aL6++N5P9e9hQ2u29/HRhzx8jf0KYUNIyrtoizYfOKZDpi/QP286WPNiU1Kma29T16K7VWXmzDm5lOzax5ozxwjac+Y0bX9C1EMCdztTVFquH3zrD7133LWnxpwP9QTjN5cm6ZDpC3RGXnGLHE+I9sSawC0DTduAv/dm8UNKKfd0jLJq/b+M/BLmfLaMnVdc32wLvjZJPRmLK5MyOCPIEx93Z9seT1K3xWlG5uNuAxIP5wKwNaA3xVFuuDRS/kBmIe+v2sN361MoKTfxbsg4vn7xLcJ/bjhwaa3JLS7Hy9WxwXInq+yJJ3GE4y7S5RWXsfHAMaYM79GsxxbidCAt7jZgy6Gc6sff//Fvvckpqe98yD3XPMHIl5fx7T8pXHFWEJ9HeeBib8cNZ95AUlpevcdIySpk8qt/MCD2Dya/upjVyRnoOXPB19cYCmcjK5OO0u/3HK69MoafXYKrVyeP351JuUlzfm8brPgtxGlOli5rA0a8vJwIpzLSdh0g2bkjc76PZWh4YPUK2gA5hWVc8cCnpLl5c+P2P7n1+9n4eRpt870ZBVw3Jx6T1nwzZQi9unhUb2cyaT5fu58XF+3ArqiIqzYvZuGZF5Du6E6/jH08tmQOQwtSISPj5E4iPp6DM19j7Fn/h7enKxrYn1mIl6sjE84JInnjDrYfK2NdeC5O/51ycscS4hRkzThuaXG3spzCMvZnFtJ/QA/uu2c8eQ4uXD/peV699mHMv1TfXbGbfV5+fPbdkzx6cFV10Abo7tuBr6cMARST3l/HnqP5gBHQJ85dS8z87QwK7cTi0Z2Jtd/Hqiu78dJVZ5Lv35XrJz1P2P99SP/Hf+f8pxfy1zVTmtRfnv7cy0zxPZ+K4mI+PsuR5Utf5KuhHpzf25cv1u5ndZELV235E6fHZmOL5sgAAAwUSURBVJz0ZybE6U76uFvZtlSjm6R/kBfn9fLlmfH9+HZ9Cm8m5VL4WyKPXxbO0fwSPlm9m/Hp2xnU1bPOSfd7dnbn69sjmTh3LZPeX8vEc4N5b8VunB3sePnqM7l6YFeUUnDhIpyBa4FxxSl89fhbZOBEUa/erHHvxv+6XcJntz1A6NuvYBo4EDSYNGg0Jm204Dt2cMLdbCHb/JVxTOwxnsNOnry7/Qe6r9oN69YxFBi6aBGZ+SUsm/M9l+z4A2bObJkPVohTmHSVtLL3VuzmhYU72PzkxXi7OUF8PDo2ltgJD/HJnhJuiAzG0d6Oz9fs4c/376D74P7HdaHUtuNILpPeWc2xUs1F/o48d+vw41rnx6laDNbHB379lfRiE1f+sItDHr4N1rlDRSmvDunE6KtGkPnXGqa+v5LVQf34YuFLDN22xkj59vY2Lk7KSA8hLCKrvLcjWw/m0K2TqxG0AWJjUYsXEwO43j+bd//aDcB1nkVG0G4knbqvvyfzNn/Kvu27GdnbF3X/JfUXrtpXZYDtAnz/VCxL8xwhNBQ1bRpKgZ1SKIx7Zr3Glw7duOMfJyapLSxZk0puQDjPrvmMoffcBPPcJWAL0cwkcLeyrYdyODPIu+aFymCqYmKYNqQPHb7+gi879OKe7b812NI21/3R++luyZwZVeOfzQS+EMuNsbHw8GQYEnLiNvdcy+VPP8vjkf35+u8Uwkty+WL9V/SlAPr3hyly4VGI5iZdJa0ou7CUs55ewvTRfblzZM+6C1V2nag21orVWrN1ws30+fVrnL29IDPTWG3Gwi8XIcTxpKukndh2yEi86R/kVX+hqChUGwyGSinOnHYnFKXDhAnG7HoyK54QLUICdyvacigbaCRwt2XmXS3SRSJEiznlxnEfyCysztZr67YdyiG4kxtebs2bgi6EOLVYFLiVUqOVUjuVUruUUo80vkXT7Fuymv/e8QbHVqxp0vap2UVc+PIyLr7vM1bO/tTGtavFfN3EJq6huOVgDv27ttPWthCi1ViyAo498DZwMXAQ+EcpNV9rnWDryhS+/R5Le19D7JdrmT1i2Anva6357Ku/+HfFRiJHDWRYZwe6vvxs9fCzHzccpAyFXUU5Nx0J5MpvN/P4ZeF02rKBAy/MxhQdTcCCebg8+fiJF/ri4yE2lpwZT7D0SBmXfvwKrk8+VucFQZNJs+7VDxi4ZClOQJa9C4uPQPzcleR/tgHf/n3w7dENX3dnfA/vx/fbz+l8+/+RW26i+OPPcLj5JvLLNQePFTE5wNafohDiVNfoqBKlVBTwlNY6uvL5owBa6+fr26bJo0ri45k1ZyGv+0fywZ5fuWjqzdWBs8KkiZm/jS/WHsCjpIA8Z2Mx2OcXvckkfzD9vpCRr/xFUOExPn73Lt6Z+hrv5nnRwQ46HT3EHu9AANxLCrl383xuf+5O1NCh1YcuuXQMP6VW8MpFt5Hh6Eb/w8l8kLmKTr98T+ryOLq+8hz2MU+ihwwh9tcEPonb9//t3XtwFfUVwPHvASQYGggEVCC8nEIt2qE4GINjVQamhShBKiAyoMwgEjp0pthx0g5Dp5S2NhVJEZ3SWEBltEhtq2mtgPIwVggvofKoFAhPGwgJkPKSEHL6x24gDY+7mL27bO75zNyZ3dxfcs+5u/dk7293fz9Gl2/hJ33bMnh7C/Ykp3HTqWO0P1FBedtbqGiRQnVN7Ct23ti8kHuWLLr298oY06j4OlmwiAwHBqnqk+76WOBuVZ1cr50vs7xXVdeQ/cOFHKI5L77/Ive0Ok/VzFn8YE8SS7YdYuJXW5D71nPsmjKVsWtOkFG6gznjMim+uSejCorJf7Q3w/qkA7Dz8AlmTF/IF8f/y5ADn9Dy4SG8u3Efy7v0IefwRnLnT0NE2LjvKJMLPqL0fDP6JFXxcM9U8j45BlpD1Q1JVEsTMg5s5dWlzzN34Dhmd7ufHilN2XniPBnH97EutSu/X5bPgE3LEfcuxJq7Mzl+5hzlQx6hfPd+jqSkcePZ06RWn+Hc7DmcnVvA0QOlfLflKZrYxLLGJDy/C/cI4Dv1CneGqn7/Sr/T0Ou493/wD0a+s4dDLdvSqbKMJK2mJLUj0x7qxfh7u19o9+Sr69lXcZplU+5jwuz3Wfv5SdZl38KN37p4JF3bBVLbnVKzejXT5hfxertvMDazK13aJpO35DM6nqrgmWUv82DrKpoUF7O9Y08WfDOL9mcqSXpiLL/51ynSTlVS/pU2jCzdxM/LVjM67QE2pN/OiIrtPHdX64uXxNXtXql9/d69Yd48Z6yOp566JC5jTGLzu3AH11VSx+m5L/Pugr/yTP+JALywdwnZi+b8X5uZS3fw2w93M/7e7hQUlfCjlQvIaXMy5k0gqsqz731GQVEJAIM73EDe4l/QanWRM85GcTHk5jqDOU2ZAnl5/GXxKuas2Mld+7bwy2PraZo/i8O503ilUwY5OQ/S+jJ98sYY45XfN+CsB3qISHfgc2AUMLoB8XmS/PafGLFuKelJSnLZIXo/PeGSNrd1SOF8jVJQVMJjXZOY2Oakp5tARIQfD76NLq/Po9nWLTzaQZCZv7p4BAyQl+c8XMM6JzFs8zxnJT8f+vXj5qqT5C56FvascIq9McYEIGbhVtVqEZkMLAWaAvNVdVvcI3MLaL/jR2DHBvhz2iU3edze0bmUrv/X2jPj8b7IpIGe/7yIMGbyIzB968XuiqsdqU+fDmvXOrd1W9eGMSZMXmcVvpaHr7O8x5jBe03hh3pmUFb8Z/i+XBw2u7gxxidcwyzv0R9kqnZMaRvgyBgTYYk1yFTdMaWNMSYBRL9wx+qbNsaYRqbRDTJljDGNXVz6uEXkCPBlbp1sB5T7HM71znJODJZzYmhIzl1Vtb2XhnEp3F+WiGzw2jnfWFjOicFyTgxB5WxdJcYYEzFWuI0xJmKut8JdEHYAIbCcE4PlnBgCyfm66uM2xhgT2/V2xG2MMSaGUAp3rDksRSRJRN50n18rIt2Cj9JfHnJ+WkS2i8inIrJcRLqGEaefvM5VKiLDRURFJPJXIHjJWURGutt6m4i8EXSMfvOwb3cRkZUissndv7PCiNMvIjJfRMpEZOsVnhcRecF9Pz4VkTt9D8LroCZ+PXBGGNwN3Ao0B/4J9KrX5nvAXHd5FPBm0HGGkHN/INldnpQIObvtUoAioBjoG3bcAWznHsAmoI27flPYcQeQcwEwyV3uBewNO+4G5nwfcCew9QrPZwHvAQJkAmv9jiGMI+4MYJeqlqhqFbAIGFqvzVCgdpr2t4ABIiIBxui3mDmr6kpVPe2uFgPpAcfoNy/bGWAG8GvgiyCDixMvOU8AXlLVYwCqWhZwjH7zkrMCrdzl1sB/AozPd6paBBy9SpOhwGvqKAZSRcTXacHDKNydgAN11g+6P7tsG1WtBiqBtECiiw8vOdc1Huc/dpTFzFlE+gCdVfVvQQYWR162c0+gp4h8LCLFIjIosOjiw0vOPwXGiMhB4O/AFac9bCSu9fN+zcIYZOpyR871L23x0iZKPOcjImOAvsD9cY0o/q6as4g0AfKBcUEFFAAv27kZTnfJAzjfqj4SkTtU9XicY4sXLzk/Bryiqs+7UyEudHOuiX94oYh7/QrjiPsg0LnOejqXfnW60EZEmuF8vbraV5PrnZecEZGBwFQgW1XPBhRbvMTKOQW4A1glIntx+gILI36C0uu+/Y6qnlPVPcAOnEIeVV5yHg8sBlDVNUALnDE9GitPn/eGCKNwX5jDUkSa45x8LKzXphB4wl0eDqxQt9c/omLm7HYb/A6naEe93xNi5KyqlaraTlW7qWo3nH79bFUNaAaOuPCyb7+NcyIaEWmH03VSEmiU/vKS835gAICIfB2ncB8JNMpgFQKPu1eXZAKVqlrq6yuEdFY2C/g3ztnoqe7PfobzwQVnw/4R2AWsA24N+0xyADl/ABwGNruPwrBjjnfO9dquIuJXlXjczgLMArYDW4BRYcccQM69gI9xrjjZDHw77JgbmO8fgFLgHM7R9XggB8ips41fct+PLfHYr+3OSWOMiRi7c9IYYyLGCrcxxkSMFW5jjIkYK9zGGBMxVriNMSZirHAbY0zEWOE2xpiIscJtjDER8z//Ma0J8Tlo2QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def lwlr(testPoint,xArr, yArr, k = 1.0):\n",
    "    xMat = mat(xArr)\n",
    "    yMat = mat(yArr).T\n",
    "    m = shape(xMat)[0]\n",
    "    #设置权重初始对角矩阵，为每一个样本点初始化了权重值为1\n",
    "    weights = mat(eye(m))\n",
    "    #遍历数据集，计算每个样本点对应权值\n",
    "    for j in range(m):\n",
    "        #样本点与待遇测点之间的距离\n",
    "        diffMat = testPoint - xMat[j,:]\n",
    "        #计算对应权值\n",
    "        weights[j,j] = exp(diffMat*diffMat.T/(-2.0*k**2))\n",
    "    xTx = xMat.T*(weights*xMat)\n",
    "    #判断是否可逆\n",
    "    if linalg.det(xTx) == 0:\n",
    "        print('cannot inverse')\n",
    "        return\n",
    "    ws = xTx.I*(xMat.T*(weights*yMat))\n",
    "    return testPoint*ws\n",
    "\n",
    "#测试函数\n",
    "def lwlrTest(testArr, xArr, yArr, k = 1.0):\n",
    "    m = shape(testArr)[0]\n",
    "    yHat = zeros(m)\n",
    "    #依次预测测试数据的yHat的值\n",
    "    for i in range(m):\n",
    "        yHat[i] = lwlr(testArr[i], xArr, yArr, k)\n",
    "    return yHat\n",
    "\n",
    "#测试一下\n",
    "xArr, yArr = loadDataSet('ex0.txt')\n",
    "#得到所有点的估计\n",
    "yHat0 = lwlrTest(xArr, xArr, yArr, 0.1)\n",
    "yHat1 = lwlrTest(xArr, xArr, yArr, 0.01)\n",
    "yHat = lwlrTest(xArr, xArr, yArr, 0.003)\n",
    "#绘制图像\n",
    "xMat = mat(xArr)\n",
    "srtInd = xMat[:,1].argsort(0)\n",
    "xSort = xMat[srtInd][:,0,:]\n",
    "fig = plt.figure()\n",
    "\n",
    "ax = fig.add_subplot(311)\n",
    "ax.plot(xSort[:,1], yHat0[srtInd])\n",
    "ax.scatter(xMat[:,1].flatten().A[0], mat(yArr).T.flatten().A[0], s = 2, c = 'red')\n",
    "\n",
    "ax = fig.add_subplot(312)\n",
    "ax.plot(xSort[:,1], yHat1[srtInd])\n",
    "ax.scatter(xMat[:,1].flatten().A[0], mat(yArr).T.flatten().A[0], s = 2, c = 'red')\n",
    "\n",
    "ax = fig.add_subplot(313)\n",
    "ax.plot(xSort[:,1], yHat[srtInd])\n",
    "ax.scatter(xMat[:,1].flatten().A[0], mat(yArr).T.flatten().A[0], s = 2, c = 'red')\n",
    "plt.show()\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "以上自下：k=0.1， k= 0.01, k = 0.003.  \n",
    "由上图知，局部加权线性回归可以更好的拟合数据。且k值越小（越精细），拟合程度越好。但是k值过小会导致过拟合。 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-12-12T00:55:07.078034Z",
     "start_time": "2018-12-12T00:55:02.993252Z"
    },
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "k=0.1, error: 56.78868743050092\n",
      "k=1, error: 429.89056187038\n",
      "k=10, error: 549.1181708827924\n"
     ]
    }
   ],
   "source": [
    "#预测鲍鱼年龄\n",
    "\n",
    "#计算平方误差\n",
    "def ressError(yArr, yHatArr):\n",
    "    return ((yArr-yHatArr)**2).sum()\n",
    "\n",
    "#不同的核大小对应的误差\n",
    "abX, abY = loadDataSet('abalone.txt')\n",
    "yHat01 = lwlrTest(abX[0:99], abX[0:99], abY[0:99], 0.1)\n",
    "yHat1 = lwlrTest(abX[0:99], abX[0:99], abY[0:99], 1)\n",
    "yHat10 = lwlrTest(abX[0:99], abX[0:99], abY[0:99], 10)\n",
    "print('k=0.1, error:', ressError(abY[0:99],yHat01.T))\n",
    "print('k=1, error:', ressError(abY[0:99],yHat1.T))\n",
    "print('k=10, error:', ressError(abY[0:99],yHat10.T))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "由上述结果显示，核越小，误差越小。但是核太小可能会导致__过拟合__，对新数据的效果差。   \n",
    "假设我们使用数据集的前100个数据作为拟合点来拟合数据集的后100个数据。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-12-12T01:11:16.190175Z",
     "start_time": "2018-12-12T01:11:12.564021Z"
    },
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "k=0.1, error: 57913.51550155911\n",
      "k=1, error: 573.5261441895982\n",
      "k=10, error: 517.5711905381903\n",
      "standard LR error: 518.6363153245542\n"
     ]
    }
   ],
   "source": [
    "yHat01 = lwlrTest(abX[100:199], abX[0:99], abY[0:99], 0.1)\n",
    "yHat1 = lwlrTest(abX[100:199], abX[0:99], abY[0:99], 1)\n",
    "yHat10 = lwlrTest(abX[100:199], abX[0:99], abY[0:99], 10)\n",
    "print('k=0.1, error:', ressError(abY[100:199],yHat01.T))\n",
    "print('k=1, error:', ressError(abY[100:199],yHat1.T))\n",
    "print('k=10, error:', ressError(abY[100:199],yHat10.T))\n",
    "ws = standRegres(abX[0:99], abY[0:99])\n",
    "yHat = mat(abX[100:199])*ws\n",
    "print('standard LR error:', ressError(abY[100:199], yHat.T.A))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "由上述结果显示，对于新数据，k=10的拟合效果较好，反而k=0.1的误差最大，也就是产生了过拟合。而且标准线性回归的效果也不比lwlr差。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.缩减系数\n",
    "缩减回归系数，也就是起到特征选择的作用。大致有两种方法，岭回归主要针对特征维度大于样本数的情况，而前向逐步回归则是一种贪心算法，通过对某特征增加或减小来判断该特征的变化是否会导致整个模型的误差的变化，因此可以起到特征选择的作用。这两种方法一般用于__回归__任务中。\n",
    "### 3.1 岭回归\n",
    "当数据的特征维度太大，大于数据个数时，由于$\\mathbf{X}^T\\mathbf{X}$ 无法求逆，导致无法使用线性回归。因此需要__系数缩减__方法来__降低特征维度__。岭回归一般用于处理__特征数多于样本数__的情况。在$\\mathbf{X}^T\\mathbf{X}$ 后加一项 $\\lambda\\mathbf{I}$ 使矩阵可变成 $\\mathbf{X}^T\\mathbf{X}+\\lambda\\mathbf{I}$, 该项可逆。回归系数公式为：  \n",
    "$$\\hat{w} = (\\mathbf{X}^T\\mathbf{X}+\\lambda\\mathbf{I})^{-1}\\mathbf{X}^Ty$$   \n",
    "__主要思想__: 加一项对角矩阵使原不可逆矩阵可逆。计算出不同组回归系数，选取在测试集上误差最小的。同时对于每一个特征，其对应的系数越大，说明该特征越应该保留，系数越小越应该放弃该特征。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-12-12T02:13:05.434059Z",
     "start_time": "2018-12-12T02:13:04.626486Z"
    },
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEKCAYAAADuEgmxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl83Hd95/HX5zeHNLqv0eHb8pHEsR3bUZyQOGcTmoQjpIQQWEpTClkotLAsu9A9CmV7sGWBpYSjuQoBWmALBAcSyOUDEpJYjh2fiS3bsi1btiQfknXN9fvsH/OTLNuyPbIl/WakzxN+j985v/n8MvK85/s7RVUxxhhjMuH4XYAxxpjcYaFhjDEmYxYaxhhjMmahYYwxJmMWGsYYYzJmoWGMMSZjFhrGGGMyZqFhjDEmYxYaxhhjMhb0u4DRVlVVpbNmzfK7DGOMySnr16/vUNXo+ZabcKExa9YsGhsb/S7DGGNyiojszWQ52z1ljDEmYxYaxhhjMmahYYwxJmMWGsYYYzJmoWGMMSZjFhrGGGMyZqFhjDEmYxPuOo0LlYileO03GZ2mPLZkjFcpp76BDPN+6Wky+OKhy4g3IiKI442LNy4gjtf3pgcCghNwcAKS7oLp4VOnOwRDDqH8AHmRIIGQM/g+xpjsYqHhScZTND7d7G8R9rh2ABxHCEUChPOD6S4SIBwZGA5SUplPWU0BZdUFlEYjBELWYDZmvFhoeCLFYT7+7Vv8LmNMqZ6WSjrMoOqQ4ZMzFD05rN66vP7AuLqnTnddRV0llVTclOKm3MF+KjUwTXGTLsmES6I/SawvSbw/RcLrx/uTxPuT9HbGOX6ol1hvkv6exGDdIlBcmZ8OEC9IymoilFUXUFSRj+NYi8WY0WShMYmcsctHhhvM/i/ZWF+SzrZejh/2urY+jh/upXVXK4lYanC5SHGIS66pY8F1dZTXFvpYsTETh4WGyTl5kSDVM0uonllyynRVpbcrTmdbL8cO9bJv61E2Pb+fjc/uo7a+lAUr6pizrJpwvv3ZG3Oh5IxdFjmuoaFB7YaFZkBvV5w3Xm5l+4utHD/cSygvwLyGai67bgo1s0vsgLsxHhFZr6oN513OQsNMBqrKoV2dbHuplabGwyTjLhVTCrns2jouuaaWSFHY7xKN8ZWFhjFnEe9LsrPxMNtfauXwni6coHDDe+dz+fVT/S7NGN9kGhq2c9dMOuFIkMuvn8rl10/lyIFuXvpZE6t/+CbxvhRL3zrD7/KMyWp2gruZ1CqnFnHnxxYz98pqXvpZE6+s3H3mqcnGmEHW0jCTXiDocNufXU4oL0DjU80k+lNc9565dpDcmGH41tIQkekiskpEtovIVhH55DDLiIj8k4g0icgmEVnmR61m4nMc4eYPXMrim6fx+gv7Wf2DN3Bda3EYczo/WxpJ4D+r6msiUgysF5FnVXXbkGXuAOZ53dXAt72+MaNOHGHFvfMIR4I0PtVMPJbi1j9dQCBge3GNGeBbaKhqK9DqDZ8Qke3AVGBoaNwFPK7pncwvi0iZiNR5rzVm1IkIV7+znlBegN//fBfJWIo/fGAhwVDA79KMyQpZ8RNKRGYBS4FXTps1Fdg/ZLzFm2bMmFr2hzO58X3zad5yhF8+uIl4f9LvkozJCr6HhogUAT8FPqWqXafPHuYlZ+xoFpEHRKRRRBrb29vHokwzCS28cRq33r+AgzuPs/LrG0+5UaIxk5WvoSEiIdKB8UNV/dkwi7QA04eMTwMOnr6Qqj6kqg2q2hCNRsemWDMpXXJ1Lbd/ZCHt+07wxNc20NsV97skY3zl59lTAjwKbFfVr55lsZXAB72zqK4BOu14hhlv9UujvO3PF9N5uJdfPvg6amdVmUnMz5bGdcAfA7eIyEavu1NEPioiH/WWeQrYDTQBDwN/7lOtZpKbcXklN33gUtr3naDptTa/yzHGN36ePfU7zvPwBu+sqY+PT0XGnNu8q2pY/3Qz637VzJxl1faAJzMp+X4g3Jhc4TjCVW+fzbHWHnatt9aGmZwsNIwZgbnLqqmYUsi6X+2xK8bNpGShYcwIiCNc9bbZHDvUS1PjYb/LMWbcWWgYM0JzlkapnFrIul8146Zcv8sxZlxZaBgzQgOtjeOHe9m5zlobZnKx0DDmAtQviVI5rchaG2bSsdAw5gKIIyx/+2w62/vY8aq1NszkYaFhzAWafUUVVdOLWPeUtTbM5GGhYcwFEkm3Nrra+3jzlUN+l2PMuLDQMOYizFpcRXRGMY1PNZOy1oaZBCw0jLkIg62Njn7e/L21NszEZ6FhzEWauaiS6pleayNprQ0zsVloGHORRITl76jnxNF+3vi93bnfTGwWGsaMghmXV1Azu4TGp621YSY2Cw1jRsHAsY3uozG2v2StDTNxWWgYM0qmL6igtr6E9U83k0pYa8NMTBYaxoySdGujnu5jMba9eMaj7I2ZECw0jBlF0y4rp25uKeufbiaZSPldjjGjzkLDmFEkkn66X09nnJ3r7Ol+ZuKx0DBmlE27pJziinx2bbDQMBOPhYYxo0xEqF8WZf+2o8T6kn6XY8yostAwZgzMWVqNm1KaN3X4XYoxo8rX0BCRx0SkTUS2nGX+TSLSKSIbve6vx7tGYy5E7ewSCkvD7N7Q7ncpxowqv1sa3wVuP88yv1XVJV73xXGoyZiLJo5QvyTKvq1HSMTsLCozcfgaGqq6FjjqZw3GjJU5y6pJJlz2bjnidynGjBq/WxqZeIuIvC4iT4vI5X4XY0ym6uaVESkO2VlUZkLJ9tB4DZipqlcA3wCeGG4hEXlARBpFpLG93fYhm+zgOMLsK6Ls3XzELvQzE0ZWh4aqdqlqtzf8FBASkaphlntIVRtUtSEajY57ncaczZylURKxFPu32V5YMzFkdWiISK2IiDe8nHS9toPY5Iypl5STVxBkl51FZSaIoJ9vLiL/BtwEVIlIC/B5IASgqt8B7gE+JiJJoA+4T1XVp3KNGbFA0GH24ip2v95BKukSCGb17zRjzsvX0FDV951n/oPAg+NUjjFjon5ZNW+8fIiWN48x8/JKv8sx5qLYzx5jxtj0y8oJ5QXY/ZqdRWVyn4WGMWMsGAowa1Elu1/vwE3Zw5lMbrPQMGYc1C+tpr87wcGmTr9LMeaiWGgYMw5mLqwkGHLYZbuoTI6z0DBmHITyAsxYWMnuje2oaycAmtxloWHMOJmzNEpvZ5xDu20XlcldFhrGjJNZi6pwgmIX+pmcZqFhzDgJR4LMuKyCXRvasGtUTa6y0DBmHNUvrab7aIy2vSf8LsWYC2KhYcw4mn1FFY4j7LbbpZscZaFhzDjKLwwx9ZIydr3WbruoTE6y0DBmnNUvraazvY8jB3r8LsWYEbPQMGac1S+JIoJd6GdykoWGMeOsoCRM3dwyO/XW5CQLDWN8MGdZlGOtPRw7ZLuoTG6x0DDGB/VLqgHY9Zq1NkxusdAwxgdF5XnUzC5hl516a3KMhYYxPpmzrJqO/d10tvf5XYoxGbPQMMYnc5ZGAay1YXKKhYYxPimpihCdUcxuO4vK5BALDWN8VL80yuE9XXQf6/e7FGMyYqFhjI9O7qKy1obJDb6Ghog8JiJtIrLlLPNFRP5JRJpEZJOILBvvGo0ZS+W1hZTXFdouKpMz/G5pfBe4/Rzz7wDmed0DwLfHoSZjxtWcpVFam47T2xX3uxRjzsvX0FDVtcDRcyxyF/C4pr0MlIlI3fhUZ8z4mLMsiirsed1aGyb7+d3SOJ+pwP4h4y3etFOIyAMi0igije3t9g/P5JbKqUWURCN2XMPkhGwPDRlm2hkPIVDVh1S1QVUbotHoOJRlzOgREeYsjXLgjWP09yT8LseYc8r20GgBpg8ZnwYc9KkWY8ZM/dIorqs0b+7wuxRjzinbQ2Ml8EHvLKprgE5VbfW7KGNGW83MEorK8+wGhibrBf18cxH5N+AmoEpEWoDPAyEAVf0O8BRwJ9AE9AJ/6k+lxowtcYT6JVG2/vYg8f4k4Xxf/2kac1a+/mWq6vvOM1+Bj49TOcb4as6yKJtWtbB3yxHmNdT4XY4xw8r23VPGTBq1c8qIFIfsQj+T1Sw0jMkSjiPMXhKlecsRkvGU3+UYMywLDWOyyJylUZKxFPu2neuaV2P8Y6FhTBaZekk5eQVBdm+0XVQmO1loGJNFAgGHWYuraN7UQSrp+l2OMWew0DAmy8xZGiXWm+TAm8f8LsWYM1hoGJNlpi+oIJQXYJftojJZyELDmCwTDAWYuaiSPRvbcd0zbrVmjK8sNIzJQvVLovSdSNDadNzvUow5hYWGMVlo5sJKAiHHbpduso6FhjFZKJwfZMaCCnZvaEdtF5XJIhYaxmSpOUuj9ByPcXhvl9+lGDPIQsOYLDVzURWOI+y226WbLGKhYUyWyi8MMe3ScnZtaCN9w2dj/GehYUwWq18apaujnyMHuv0uxRggw9AQke+LyEdE5NKxLsgYc9LsK6KIYE/0M1kj05bGvwB1wDdEZJeI/FREPjmGdRljgIKSMHVzy+zUW5M1MgoNVX0B+DvgfwKPAA3Ax8awLmOMZ86yKMdaezh2qMfvUozJePfU88CLwHuBN4GrVNV2VRkzDuqXVANYa8NkhUx3T20C4sBCYDGwUEQiY1aVMWZQUXkeNbNL7DGwJitkunvqP6nqDcDdwBHSxzjspjjGjJP6pVHa952gq6PP71LMJJfp7qlPiMiPgY3Au4DHgDsu9s1F5HYReVNEmkTkc8PMv19E2kVko9d9+GLf05hcNHdZNQhsWXvA71LMJBfMcLkI8FVgvaomR+ONRSQAfBO4DWgB1onISlXddtqiP1bVT4zGexqTq0qqIsxrqGHz6haW3jaDSHHY75LMJJXp7qkvq+oroxUYnuVAk6ruVtU48CPgrlFcvzETylVvm0Uq4bLh2X1+l2ImMT+vCJ8K7B8y3uJNO927RWSTiPy7iEwfn9KMyT7ltYXMuyrd2ujtivtdjpmk/AwNGWba6TfYeRKYpaqLgeeA7w27IpEHRKRRRBrb2+0MEzNxNdyZbm1stNaG8YmfodECDG05TAMODl1AVY+oaswbfRi4crgVqepDqtqgqg3RaHRMijUmGwy2NtZYa8P4w8/QWAfME5HZIhIG7gNWDl1AROqGjL4T2D6O9RmTlQZaG3Zsw/jBt9DwDqp/AvgN6TD4iapuFZEvisg7vcX+UkS2isjrwF8C9/tTrTHZo7y2kHnLa9hixzaMD2Si3ae/oaFBGxsb/S7DmDF1/HAv//qFl7ni1hlc9+65fpdjJgARWa+qDedbzp6nYUwOKqspYP7yWmttmHFnoWFMjmq4cxappMuGZ/b6XYqZRCw0jMlRg62NNQestWHGTaa3ETHGZKGGO2ex49VDrHu6ifk3lHP06FGOHj1KPB4nLy+P/Px88vLyhh0Oh8M4jv1uNCNjoWFMjujt7eXIkSMcO3ZsMByOHj3Ksbp2Vm2PsWrICekiQiYnueTn5zN9+nTmzp3L3LlzqaysHMMtMBOBhYYxWSqVStHS0kJTUxM7d+7k0KFDp8wvKSmhoqKCSy65hD2vnqB+wTSuvXMBFRUVhMNhEokE/f39xGIxYrHYsMPd3d3s2bOHnTt3AlBeXs68efOYO3cus2bNIhy2GyOaU1loGJNFTpw4MRgSu3fvpr+/HxFh+vTp3HLLLdTU1FBRUUFZWRmhUGjwdc/1bWPX+jZK7qkgLy8PgHA4nPGX/tGjRwff97XXXuPVV18lEAgwc+ZM5s6dy7x586iqqkJkuLv/mMnErtMwxkcDrYmdO3fS1NQ02JooKioa/MVfX19PJHLuB2Ueb+vlX7/wCotvmcaKe+ZdVE2JRIJ9+/YNhkhHRwcA0WiUFStWsHDhQgKBwEW9h8k+mV6nYaFhzDjr7e2lqamJHTt20NTUdEprYiAoamtrR/yr/vnvbqNpfRsf+Nu3UFiaN2r1Hj9+nJ07d7Ju3Tra2tooLy/n+uuvZ/HixQSDtrNiorDQMCZLqCodHR3s2LGDHTt2sG/fPlSVgoIC5s+fz7x58zJqTZzPYGvj5mmseM/FtTaG47oub775JmvXrqW1tZXS0lJWrFjBkiVLTtlVZnKThYYxPkomk+zdu3cwKI4dOwZATU0N8+fPZ/78+UydOnXUT3l9/nvb2NnYxh+PcmtjKFWlqamJNWvW0NLSQnFxMddeey1XXnmlHTjPYRYaxoyTVCpFe3s7ra2tg92hQ4dIJBIEAgHq6+sHg6K0tHRMaxlobSy8YSo33Dd/TN9LVdmzZw9r166lubmZgoICrr32Wq666qrBg/Emd2QaGrZD0pgRSCQStLW1nRIQhw8fJpVKARAKhaitrWXp0qXMmTOH2bNnj+uv77LqAhasmMLm1S3kF4W46m2zxuyMJxGhvr6e+vp69u7dy9q1a3nuued48cUXectb3sLVV19t4TEBWUvDGNJh0N3dPdj19PScMj7QdXZ2Dl40l5+fT11dHXV1ddTW1lJXV0dlZaXvV1mnUi6rv/8Gb7x8iIU3TuX6987HccbnVNkDBw6wZs0aduzYQSQS4dprr2X58uUWHjnAdk+NUDKZZNeuXWNQ0cST6d/M6csNHR9u3sC08w2rKq7rnrUbmJ9MJkkmkyQSibN2yWSSeDxOPD78vZsKCgooKioa7EpLSweDoqysLGuvW1BVfv/zXWx4Zh9zlkW57U8vJxAavzA7cOAAq1evZufOnYO7rZYvX27HPLKYhcYI9fT08OUvf3kMKjLjzXGcwS4cDhMKhQgGg4RCobN2hYWFFBYWnhIQhYWFOX89woZn9/HST5uYekk5d350EeHI+O6RbmlpYdWqVezatYuCggJWrFhBQ0ODhUcWstAYof5EP6/ueHUMKsodF/2reZiXy8BEOfe0wfcXb/7QvlfXwP2UxBEcx0FEEEcQ8ca94QHKmX/bw047vdWDnjF96LTThxU9OV0ZnHb68kPXOfCaAY6kt8cRB4chw+IMznfEIRwIk+fkkRfMIz+QTzgQJj+YT9gJn/Xze/PlVl54/A0qpxXx9k9cQUHJ+H9h79+/n1WrVrF7924KCwsHw8NO1c0eFhojdLT/KDf++MYxqMiY8ZEXyEuHSCCfglABZXllg11F2wzyXqgnUKjMeJ8QrS2lPK+csvwyKvIrBsNprO3du5fVq1ezZ88eioqKuOaaa1i4cCFlZWXj8v7m7Cw0RiieivNy68tjUJHxkwzT/BnuF/nAcqe3goa+frDFM/C/oeMip/SHW35g2tD3EGSwheKqi4t3TEZdXHVRdLCfclPEU3FiboxYMkYsNUznTe9J9HA8dpzOWCfHYsfojHVScqyGO7Y/QMpJ8tRl3+FI4UEAgk6QmoIa6grrqCuso7awltrC2sHxuqI6CkOFF/MxnKG5uZnVq1fT3NwMwLRp07j88stZsGDBmJ+WbIZnoWGMOUV/sp99ew+z9qFmEv0pat+doL/6KId7DtPa08qhnkO09rTS1ttGSlOnvLY4XExdYR1TiqYwrWgaU4qmnDJcHC6+oJqOHj3K1q1b2bp16+B9t2bMmDEYIMXFF7ZeM3IWGsaYYZ042s+T/7SRro5+3vpnl1O/NHrK/KSbpKOvYzBEWntaOdh9kEM9hzjQfYAD3QfoS/ad8pricPEpYVJbUEtNYQ01BTXUFtZSFaki6Jz7IHxHR8dggLS1tQEwc+ZMFi5cyGWXXUZRUdHo/ocwp7DQMMacVX93gl9+83XamruYdlkFsxZVMmtRFSVV57//lapyPHacg90HOdB9YLA/MHyw5+AZoeKIQ1V+1WCQ1BTWUFtQS1VBFRX5FVTmV1KRX0FZfhkhJ0RbW9tggAzcZbeiooLq6mqqq6upqamhurqaioqKnD/DLVvkRGiIyO3A14EA8Iiqfum0+XnA48CVwBHgvarafK51XmhodPYl+Oj314/4dWb8XfRJXuc4y+v0eUOPf8jgNHAk/Yr0bME78cvrp8cdR3BECHjLp8ch4E13RAaHQ0Eh5DgEA0Io4BAKCEHHIRR0CDlC0JsWCQUoCAeJhAMUeF16OEgkFCAwgov4ErEUjU81s3tjO8cP9wJQMaWQWYurmL24iupZJRd0UaCq0tl/nMPHW2jvPEh7VytHug5xtKuNY91tHD/RTmd3B6lYP8EU6c5lcLhYIpQ4BZRIhEKJEHJK6Jcy+p18eiVML4HBD8pBKQMqESqAClcpFaEgGEw/zjYcRkKhdHf6cDiEU1hIoKQEp6iIQEkJgeJipKAga6+/GUtZfxsREQkA3wRuA1qAdSKyUlW3DVnsz4BjqjpXRO4D/jfw3rGoR2MxFm3+7Vis2oyi4U6ZHYWVnrJuPW360EHR9FLqnVo7+Fo9/bWKq4rqwPDJackhwwqom74YMeWePA13MKBQZEhB4tVwcp56y6b7QYFwQMgPCgVBh0jQIRIUIkEhP+AQCUBeUIgEHMIBqHSEKY4Sq8yjtb+cQ22VbPj1CV779V7CxKiRQ9S4LURTBwkk+9FEIt3F4yeHT+/icUgmAaj0upHp8bpTpQRcB5KBAJ0lxXSWldFVWsqJklL2lZSyM1IAAyGXSuJ0x8jr7yc/FiO/v5+8gX5/LD0c6yeYSBJKJAgmkwST6WFHhMBAmBQXEyguxikpJlBaSqCkNN0vLSFQWoozMK2sdHB5meDPXffz3lPLgSZV3Q0gIj8C7gKGhsZdwBe84X8HHhQR0TFoHhW5ce5+9rHRXq0xWcFFUBFcBFfSwwlxiCOkxCHlOASdAFMDAWrChZwou4Su0vkcLJrD/sBMJJCk0O0gGI4RDMUIFcQIaYwQfYS0j7D2EtZuwtpNKNWFkEAcwOvS19QMjHvDcnIcSbfO1En38WoUJ4hIMUoJrpOPSh5KmELJI0KYWskDwghhkokAPbj0iUuMJDFJEstPEo+k6CXJcUkQI4kr5/76EIWgOgRUCA50ruCkIHAkQaC9A8dtx1El4CpOysVRr++6SPqngHd2nNcXHdxucRwION6wIIGA99/CQQLO4HzHcSAQwAk44Dio47VivWNDEnDAa6mm/xs6REpLuPbe94zZ3xH4GxpTgf1DxluAq8+2jKomRaST9A+XjtEuJlBaypznnhvt1ZocMaK9EcPu3zrPNDnZdjg5behsObnM2fpDlxvacfKLdmgnTvrLhlQMifeQinXTc+IY3V3H6e3upL+7i1hvJ9JzlGDvUaT/GKHYcfISx5mVXENR6kmK3F6OJOayJ3YVx5LTSWiEuEbodsuIa4S4FqJc+C/rfIWIChGBiJzsFziQ7wj5p+8e05N9VSWhLklNknQTJDROvttPyo2jBME7VRncIS9zSeISF5e445JyNN2CEUg5kJL0eEqUpCgpcdP9QLofQ3FFSeGSEsXFJYmgAum97BdzfMVNdwokL2wN5XsjXHvvRZSQAT9DY7h/pqf/BMhkGUTkAeABSJ+ud0HFJLoJr/nUBb3WjLcxOKjB6V/QZ35hn/p6GWGf06Zx6jz1fp0qXv+Ub8f0tRvxfFLxAJoI4CYCuMkgbiKIJkO4yTBuKoybysN188EVRPsR7UOIIcQRiQNxCohTSAyRBEIMR7pxnF6cPMEpCOKUVOMUz0QKyqCggpq8coqDJfQGSukPldIbLKPHKaE3UExvUujvT9LfmyTWl6C/L0miP4UqBFJKQcIlknCJxF0KEl4XT3eRxJlxk3CgLy9AX9jheJ5Db9Clr7eNnq59xPqOkOg9SqKng0Syl5QmTr4wGMYpLMcpLMMpKkaCYQiEkEAICYaQQBiC3vjANCcIzskveUW8BsLQz1wIKgRV000Qd2B/pJseVkVcRZJJxE0gqQSoi7gp0CTquqi6iNdHXe8jTbdIVN307kevr3i7H4d8/EP2nw7+ZSpeK+a02cGge+bf6yjzMzRagOlDxqcBB8+yTIuIBIFS4OjpK1LVh4CHIH0g/IKqURd62i/opWYcXfSeyWFer3rqvGEPapz6JT6i/uCqzjZPB8MkpUUkUzUkU7XpfrKaZKqaZLIa5ey3/xAnhhOI4wTiSDiBBEC1DFdDqAbTneugKQdN/5Q+cyUxoGvIOvMC6RApCOFEggRDDkUKha4S1Ti4Hair6R/IqulhBU26uCfiuL2n/Vx2hEBpmEB1hGBZPoGyvMEuWJruS16AIy372LOhkQMb1nHwze24qRThSAEVU6ZSNb2W4qpFlESrKa6KUlJVTUlVlPyi4kl58NoPfobGOmCeiMwGDgD3Ae8/bZmVwJ8AvwfuAV4Yi+MZAETK4YHVY7JqY4bSlEvySD/Jtl4S7b0k2/tIdqS7U75oHSFYkU+wKkJeVYRgVT5OYRgnEsDJD+JE0p3kBZHAyL4wVRWSiiZSuH3JdNebxO1NnDrcOzAvgduTSO87HzjuMHAsIgg4TvpMK2//emB2aToMyvIIlHsBURxGhjkbKxHrZ//Wzez+5Tr2bGykqz19jUZ0xiwa3vFHzF7awJR5l+LYqbVZwbfQ8I5RfAL4DekdgY+p6lYR+SLQqKorgUeB74tIE+kWxn1+1WvMSLmxFMn2XhLtfemAaOsl2dZL8kj/yd0cQKAkTDAaIbKoimBVAcFohGBVhGB5Xvpg5xgQEQgJEnJwCsb/poGqyo6Xf8fW1c+xf+tmkok4wbw8Zi5awtXvupfZSxsorqwa97rM+dnFfcZcADeWItUVI9UVT3edMdyu+MlpnTFSnUOe0eFAsDJCMFpAqLqAYHWEUDTdd/Im1wM0D+54g9WPP0zrzjcprallzrLlzF7awLQFiwjaXW99k/XXaRgzWnTgwOHJCyjS01xFUwN9F1JDxpNueh98Kj1PYy5uPIXGkunhWBKNpXBjKTSeGhx2exOkOuNoLHVGHZIXIFASJlCaR96cMoJVkXRARCMEKyNIcGKfv38+Xe1trP3X7/LmS2spLCvnDz/6SRbceAuOY7udcomFhifVk+Dw1ybxFeEjanBmuPA5FjujgTt4oPjc03TotKFhMdoEJBxIHwzO8/rhAKFoAflzy3FKwumAKMlLH9wtCU+6FkOm4n29vPLE/2P9r55AEK55931c9c53E84//y1LTPaxv3KPBIXI5SO/dnVCGYuzT86xymHPdjmkyG97AAAP1klEQVTlbMdTb1M+cJrqyWsSGLx3x+A057TxgKQvmAp4B3ADJ6cR8B7iFJAhARFE8gJIyBn2oK3JnOum2LLqWV788Q/o7TzOZdffzIr7PkhJVfT8LzZZy0LD4+QFKb97nt9lGDMh7N20kdXff4SOfc1MuWQB7/qv/5O6uZf4XZYZBRYaxphR03P8GM889A12r3+VkmgNb//U55h/zXV2DcUEYqFhjBkVrU1vsvIrf09/dzfXv/9+lt3xToLh8X8euRlbFhrGmIu2ZdWzPPfotygsq+B9/+vLVM+q97skM0YsNIwxFyyVTLL68UfY+JtfMmPhFbz9U58lUlzid1lmDFloGGMuSG/ncZ782pdo2b6FK99+Nze8/3671cckYKFhjBmxQ7t28ouv/B39J05w5198hstW3OR3SWacWGgYY0Zk65rnefbhByksK+e+L/4jNbPn+F2SGUcWGsaYjKSSSdb84FE2PP0k0y9fzNs/9VkKSkr9LsuMMwsNY8x59XZ18uTX/oGWbVtYdudd3PiBD9nxi0nKQsMYc06Hmnaw8qv/QF9XJ3d8/NMsuOEWv0syPrLQMMac1eYXnuH5R79FYXlF+vhF/Vy/SzI+s9AwxpwhmUiw6l/+mU3P/5oZi5bwtr/8L3b8wgAWGsaY05w40sHKr/49h5p2sPyue7juvj+2Z16YQRYaxphB+7dt5pf/93+TiMV4x6f/ivlXX+d3SSbLWGgYY1BVXnvqF6z5wWOU1U7h3r/+ByqnTfe7LJOFLDSMmeQS/f0889A3eOPFNcy96hpu//NPk1dQ4HdZJktZaBgziR0/1MovvvJ3dOzfy4r7Psjyu+5BnMn9LHNzbhYaxkxCqsqWVc+y5vuPIiK8+3NfYNaSK/0uy+QAX0JDRCqAHwOzgGbgXlU9NsxyKWCzN7pPVd85XjUaM1EdObCf5x7+Ji3btzD10gXc8fFPU1pd63dZJkf41dL4HPC8qn5JRD7njX92mOX6VHXJ+JZmzMSUTCR49Ymf8OoT/49gXh63PfAXLLr5NtsdZUbEr9C4C7jJG/4esJrhQ8MYMwr2b9vMsw9/k2MHW7j0uhu56YMfprCs3O+yTA7yKzRqVLUVQFVbRaT6LMvli0gjkAS+pKpPDLeQiDwAPAAwY8aMsajXmJzUd6KLNT94jK2rn6O0uoY/+qu/YbYduzAXYcxCQ0SeA4bbUfrfR7CaGap6UETqgRdEZLOq7jp9IVV9CHgIoKGhQS+oYGMmEFVl++9Ws/rxR+jvPsFVd93DW959H6G8fL9LMzluzEJDVW892zwROSwidV4row5oO8s6Dnr93SKyGlgKnBEaxpiTjh9q5blHv8XeTRuom3sJt/2PvyU6c7bfZZkJwq/dUyuBPwG+5PV/cfoCIlIO9KpqTESqgOuAfxzXKo3JIYlYP6/+4t9Zt/KnBIIh/uBDH2PxbbfbfaPMqPIrNL4E/ERE/gzYB7wHQEQagI+q6oeBy4B/FhEXcEgf09jmU73GZC1VpanxZVZ/72G62tu49LobufEDH6KootLv0swE5EtoqOoR4A+Gmd4IfNgbfglYNM6lGZNTjrUe4IXvPkTzxvVUTZ/JvZ//B6YvsH82ZuzYFeHG5KBEfz+vPPETGp/8GYFQmJs++BGW/OHbCATtn7QZW/YXZkwOUVWaXv09qx5/mBMd7Sy4/mZu+MCH7JoLM24sNIzJEUcPHuCFf/kOezdtIDpjFnf+zWeYdunlfpdlJhkLDWOynOumaHzy57z0kx8QCIW5+f7/yJK33okTsLOizPiz0DAmix07dJBff/NrHNyxnXlXX8sffOhjtivK+MpCw5gspKq8/uzTrPnBowSCQe78i89w6XU3IiJ+l2YmOQsNY7LMiSMd/OY7X2fvpg3MumIZb/3oX1JcUeV3WcYAFhrGZA1V5Y3freb5f/kOqWSSWz/85yy+9Q5rXZisYqFhTBbo7erkuUe+yc5XXmLK/Mu4/eP/ifLaKX6XZcwZLDSM8dmu9a/wzD9/g1hPN9e//34a3nG33S/KZC0LDWN80tvVydofPMbWNc8TnTmbe/7H3xKdMcvvsow5JwsNY8aZqrJt7Qus/v6jxHt7uPru9/KWe+4jEAz5XZox52WhYcw4OtZ6gOce+Rb7trzOlPmXcdtHPk6VtS5MDrHQMGYcpJIJ1q38GS//7EcEQ+H0mVF/cDviOH6XZsyIWGgYM8YOvLGNZx9+kCMt+5h/zQpuvv8Bisor/C7LmAtioWHMGOnv6ea3P/wum57/NcVVUe7+7OepX3aV32UZc1EsNIwZZa6bYsfvf8fqxx+ht7OTK9/2Lq699z8Qzo/4XZoxF81Cw5hRcuJoB1tWPcvmF57hREc7NfVzufuzn6emfq7fpRkzaiw0jLkIbirFno2NbHru1+zZsB5VlxmLlnDjBz7EvOXX2u3LzYRjoWHMBehqb2PzqmfY8sIzdB87SkFpGVfd9W4W3fxWymrr/C7PmDFjoWFMhpKJBHteW8emF35D8+uvATD7imXc8qGPUr9suT2f20wKvvyVi8h7gC8AlwHLVbXxLMvdDnwdCACPqOqXxq1IM6n1dZ+gvXkP7Xt309a8m/bm3Rw5sB83laKoopJr/ug+Ft18GyXRar9LNWZc+fXTaAvwR8A/n20BEQkA3wRuA1qAdSKyUlW3jU+JZqJTVRL9ffR0Hqdj/17am3fT1ryHtuZdnOhoH1yusLyC6pmzqb9yOVMuuYxZi5fZsQozafkSGqq6HTjfcwKWA02quttb9kfAXYCFxiSlrovrurjJJMlEnGQiTiqRJJWIk0wkSCXipBIJbzhBvL+Pvq4u+ru76DvRRV9XF33dJ9LDJ7roP9FFKpkcXL+IQ/mUqUy9ZAHRt86melY90Zmz7fGqxgyRzTthpwL7h4y3AFeP1Zv1dZ/gx5//7FitPiOq6uebnznpvMvqkElDltaBnoKqt7iml9GT01FNL6rpea7rghcM6nWu66KaHr5QIg75RUVEikuIlJRQWl1L7Zz5REpKiBQVEykppXLadKqmzySUl3/B72PMZDBmoSEizwG1w8z676r6i0xWMcy0Yb/HROQB4AGAGTNmZFzjUI7jUDl1+gW9dlT5+ZS2Yd77rNV4y56ttTg4XSS9DhFvmqT/7w2ne+l54gRwHMcbdhDHSY8PdJIeD4RCBEIhgqHw4PDgeDBEMBwiEAwRyo8QKSkhv6DQ7vFkzCgZs9BQ1VsvchUtwNBv8WnAwbO810PAQwANDQ0X9HM9r6CQd3z6ry7kpcYYM2lk88+vdcA8EZktImHgPmClzzUZY8yk5ktoiMjdItICvAX4lYj8xps+RUSeAlDVJPAJ4DfAduAnqrrVj3qNMcak+XX21M+Bnw8z/SBw55Dxp4CnxrE0Y4wx55DNu6eMMcZkGQsNY4wxGbPQMMYYkzELDWOMMRmz0DDGGJMx8fXWFWNARNqBvRexiiqgY5TKyQYTbXtg4m3TRNsemHjbNNG2B87cppmqGj3fiyZcaFwsEWlU1Qa/6xgtE217YOJt00TbHph42zTRtgcufJts95QxxpiMWWgYY4zJmIXGmR7yu4BRNtG2BybeNk207YGJt00TbXvgArfJjmkYY4zJmLU0jDHGZMxCwyMit4vImyLSJCKf87ue0SAizSKyWUQ2ikij3/WMlIg8JiJtIrJlyLQKEXlWRHZ6/Zx6FutZtukLInLA+5w2isid51pHNhGR6SKySkS2i8hWEfmkNz0nP6dzbE8uf0b5IvKqiLzubdPfeNNni8gr3mf0Y+8RFOdfn+2eAhEJADuA20g//Gkd8D5VzennkYtIM9Cgqjl5frmI3AB0A4+r6kJv2j8CR1X1S164l6uqv8/pHYGzbNMXgG5V/T9+1nYhRKQOqFPV10SkGFgPvAu4nxz8nM6xPfeSu5+RAIWq2i0iIeB3wCeBTwM/U9Ufich3gNdV9dvnW5+1NNKWA02qultV48CPgLt8rmnSU9W1wNHTJt8FfM8b/h7pf9A54yzblLNUtVVVX/OGT5B+9s1UcvRzOsf25CxN6/ZGQ16nwC3Av3vTM/6MLDTSpgL7h4y3kON/KB4FnhGR9d5z1CeCGlVthfQ/cKDa53pGyydEZJO3+yonduWcTkRmAUuBV5gAn9Np2wM5/BmJSEBENgJtwLPALuC497A7GMF3noVGmgwzbSLst7tOVZcBdwAf93aNmOzzbWAOsARoBb7ibzkjJyJFwE+BT6lql9/1XKxhtienPyNVTanqEmAa6T0rlw23WCbrstBIawGmDxmfBhz0qZZR4z0JEVVtI/2kxOX+VjQqDnv7nQf2P7f5XM9FU9XD3j9qF3iYHPucvP3kPwV+qKo/8ybn7Oc03Pbk+mc0QFWPA6uBa4AyERl4emvG33kWGmnrgHne2QRh4D5gpc81XRQRKfQO5CEihcBbgS3nflVOWAn8iTf8J8AvfKxlVAx8uXruJoc+J+8g66PAdlX96pBZOfk5nW17cvwziopImTccAW4lfaxmFXCPt1jGn5GdPeXxTqH7v0AAeExV/87nki6KiNRz8jnsQeBfc22bROTfgJtI343zMPB54AngJ8AMYB/wHlXNmQPLZ9mmm0jv9lCgGfiPA8cDsp2IrAB+C2wGXG/yfyN9HCDnPqdzbM/7yN3PaDHpA90B0g2Fn6jqF73viB8BFcAG4AOqGjvv+iw0jDHGZMp2TxljjMmYhYYxxpiMWWgYY4zJmIWGMcaYjFloGGOMyZiFhjEZEJHu8y+V0Xq+ICKfyWC574rIPedbzpjxZqFhjDEmYxYaxoyAiBSJyPMi8pr3rJK7vOmzROQNEXlERLaIyA9F5FYRedF7XsHQ205cISIveNM/4r1eRORBEdkmIr9iyA3+ROSvRWSdt96HvKuWjfGFhYYxI9MP3O3dCPJm4CtDvsTnAl8HFgOXAu8HVgCfIX1V8YDFwNuAtwB/LSJTSN+a4hJgEfAR4Nohyz+oqld5z9+IAG8fo20z5ryC51/EGDOEAH/v3THYJX076Rpv3h5V3QwgIluB51VVRWQzMGvIOn6hqn1An4isIn3zuxuAf1PVFHBQRF4YsvzNIvJfgQLSt3zYCjw5ZltozDlYaBgzMv8BiAJXqmrCezpivjdv6H173CHjLqf+Wzv93j16lumISD7wLdJPYNzvPeUv//TljBkvtnvKmJEpBdq8wLgZmHkB67jLe25zJembFa4D1gL3eQ/LqSO96wtOBkSH94wHO6PK+MpaGsaMzA+BJ0WkEdgIvHEB63gV+BXpO8D+L1U9KCI/J/34zc2kn1e/BtLPPxCRh73pzaQDxhjf2F1ujTHGZMx2TxljjMmYhYYxxpiMWWgYY4zJmIWGMcaYjFloGGOMyZiFhjHGmIxZaBhjjMmYhYYxxpiM/X+3MsBPxkCnNgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#岭回归主函数\n",
    "def ridgeRegres(xMat, yMat,lam=0.2):\n",
    "    xTx = xMat.T*xMat\n",
    "    #加一项对角矩阵，使整体可逆\n",
    "    denom = xTx + eye(shape(xMat)[1])*lam\n",
    "    #仍然需要判断是否可逆，因为当lam=0时，可能不可逆\n",
    "    if linalg.det(denom) == 0:\n",
    "        print('cannot inverse')\n",
    "        return \n",
    "    ws = denom.I*(xMat.T*yMat)\n",
    "    return ws\n",
    "\n",
    "#岭回归测试函数\n",
    "def ridgeTest(xArr,yArr):\n",
    "    xMat = mat(xArr)\n",
    "    yMat=mat(yArr).T\n",
    "    yMean = mean(yMat,0)\n",
    "    #数据标准化\n",
    "    yMat = yMat - yMean     \n",
    "    xMeans = mean(xMat,0)   \n",
    "    xVar = var(xMat,0)\n",
    "    #减均值除以方差\n",
    "    xMat = (xMat - xMeans)/xVar\n",
    "    #使用30个不同的lam值进行测试\n",
    "    numTestPts = 30\n",
    "    wMat = zeros((numTestPts,shape(xMat)[1]))\n",
    "    for i in range(numTestPts):\n",
    "        ws = ridgeRegres(xMat,yMat,exp(i-10))\n",
    "        wMat[i,:]=ws.T\n",
    "    return wMat\n",
    "\n",
    "#进行测试\n",
    "abX,abY = loadDataSet('abalone.txt')\n",
    "ridgeweights = ridgeTest(abX, abY)\n",
    "#查看缩减效果\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111)\n",
    "ax.plot(ridgeweights)\n",
    "plt.xlabel('lambda')# make axis labels\n",
    "plt.ylabel('w')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "由上图可知，随着lamda的取值不同，回归系数w也在变化。lambda极小时，回归系数和标准线性回归一样，当lambda极大时，所有的回归系数都趋于0.最佳的lambda趋于之间。为了获得最佳lambda，一般会进行__交叉验证__。  \n",
    "### 3.2 前向逐步回归\n",
    "基于贪心的回归系数缩减方法。  \n",
    "__主要思想__：对每一个特征增大或减小，观察对误差的影响。最后的结果是部分特征不论增大还是减小都对误差影响不大，则可以省去。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-12-12T02:40:15.674186Z",
     "start_time": "2018-12-12T02:40:14.142232Z"
    },
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "回归系数更新过程：\n",
      "[[ 0.04  0.    0.09  0.03  0.31 -0.64  0.    0.36]]\n",
      "[[ 0.05  0.    0.09  0.03  0.31 -0.64  0.    0.36]]\n",
      "[[ 0.04  0.    0.09  0.03  0.31 -0.64  0.    0.36]]\n",
      "[[ 0.05  0.    0.09  0.03  0.31 -0.64  0.    0.36]]\n",
      "[[ 0.04  0.    0.09  0.03  0.31 -0.64  0.    0.36]]\n",
      "[[ 0.05  0.    0.09  0.03  0.31 -0.64  0.    0.36]]\n",
      "[[ 0.04  0.    0.09  0.03  0.31 -0.64  0.    0.36]]\n",
      "[[ 0.05  0.    0.09  0.03  0.31 -0.64  0.    0.36]]\n",
      "[[ 0.04  0.    0.09  0.03  0.31 -0.64  0.    0.36]]\n"
     ]
    }
   ],
   "source": [
    "#数据标准化函数\n",
    "def regularize(xMat):\n",
    "    inMat = xMat.copy()\n",
    "    inMeans = mean(inMat,0)   \n",
    "    inVar = var(inMat,0)\n",
    "    #减均值除以方差\n",
    "    inMat = (inMat - inMeans)/inVar\n",
    "    return inMat\n",
    "\n",
    "def stageWise(xArr,yArr,eps=0.01,numIt=100):\n",
    "    xMat = mat(xArr); yMat=mat(yArr).T\n",
    "    #数据标准化\n",
    "    yMean = mean(yMat,0)\n",
    "    yMat = yMat - yMean     \n",
    "    xMat = regularize(xMat)\n",
    "    m,n=shape(xMat)\n",
    "    #用以记录所有迭代中产生的回归系数，用作分析\n",
    "    returnMat = zeros((numIt,n))\n",
    "    ws = zeros((n,1)); wsTest = ws.copy(); wsMax = ws.copy()\n",
    "    #迭代numIt次\n",
    "    for i in range(numIt):\n",
    "        if i > 190:\n",
    "            print (ws.T)\n",
    "        lowestError = inf; \n",
    "        #对每个特征\n",
    "        for j in range(n):\n",
    "            #增加（1），或减小（-1）该特征的值，eps指的是特征变化的步长。\n",
    "            for sign in [-1,1]:\n",
    "                wsTest = ws.copy()\n",
    "                wsTest[j] += eps*sign\n",
    "                yTest = xMat*wsTest\n",
    "                rssE = ressError(yMat.A,yTest.A)\n",
    "                #若改变该特征值后误差降低，则更新该特征值。\n",
    "                if rssE < lowestError:\n",
    "                    lowestError = rssE\n",
    "                    wsMax = wsTest\n",
    "        ws = wsMax.copy()\n",
    "        returnMat[i,:]=ws.T\n",
    "    return returnMat\n",
    "\n",
    "#测试上述代码\n",
    "xArr, yArr = loadDataSet('abalone.txt')\n",
    "print('回归系数更新过程：')\n",
    "returnMat = stageWise(xArr, yArr, 0.01, 200)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "上面打印了后10轮每轮回归系数更新的结果，可以看出__第一个__和__第六个__特征系数都为0，说明这两个特征对整体的回归任务不影响，可以省去。\n",
    "\n",
    "## 4.权衡方差和偏差\n",
    "__方差__：数据扰动对模型造成的影响，也指模型对不同数据的适应能力。模型越简单，方差越小，因为不论在已知还是未知数据集上表现都不会特别好。  \n",
    "__偏差__：模型拟合训练数据的能力。模型越复杂，偏差越小，在训练集上表现越好，未知测试集上可能越差（无法适应未知数据集）。   \n",
    "方差度量方法：取两组数据分别拟合回归模型，得到的两组回归系数之间的差异即表示了方差的大小。   \n",
    "参考：https://blog.csdn.net/simple_the_best/article/details/71167786   \n",
    "> 偏差-方差窘境 (bias-variance dilemma).   \n",
    "给定一个学习任务, 在训练初期, 由于训练不足, 学习器的拟合能力不够强, 偏差比较大, 也是由于拟合能力不强, 数据集的扰动也无法使学习器产生显著变化, 也就是欠拟合的情况;    \n",
    "随着训练程度的加深, 学习器的拟合能力逐渐增强, 训练数据的扰动也能够渐渐被学习器学到;   \n",
    "充分训练后, 学习器的拟合能力已非常强, 训练数据的轻微扰动都会导致学习器发生显著变化, 当训练数据自身的、非全局的特性被学习器学到了, 则将发生过拟合.\n"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [default]",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
